#Load data
rm(list=ls())
library(readxl)
mydata<-read_excel("Data_summary_NP.xlsx",sheet = "Sheet2")
head(mydata)
## # A tibble: 6 × 33
##   Forest     N     P Treatment Total_C Total_N TC_TN  WSOC   NH4   NO3 availP
##    <dbl> <dbl> <dbl> <chr>       <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1      1     1     1 CK           33.3    2.24  14.9  55.6  3.53  13.4  0.309
## 2      1     1     1 CK           29.2    1.72  17.0  49.9  2.80  10.2  0.275
## 3      1     1     1 CK           34.8    2.22  15.7  62.6  1.78  16.6  0.710
## 4      1     1     1 CK           37.2    2.21  16.8 107.   1.97  15    0.615
## 5      1     1     1 CK           29.3    1.83  16.0  46.0  2.99  14.2  0.553
## 6      1     2     1 N            22.9    1.61  14.2  39.6  1.39  11.0  1.11 
## # ℹ 22 more variables: availS <dbl>, pH <dbl>, EC <dbl>, MBC <dbl>, MBN <dbl>,
## #   MBC_MBN <dbl>, ARS <dbl>, SARSmbc <dbl>, SARSt <dbl>, BG <dbl>, SBG <dbl>,
## #   NAG <dbl>, SNAG <dbl>, XYL <dbl>, SXYL <dbl>, LAP <dbl>, SLAP <dbl>,
## #   ACP <dbl>, SACP <dbl>, logS_C <dbl>, logS_N <dbl>, logS_P <dbl>
str(mydata)
## tibble [60 × 33] (S3: tbl_df/tbl/data.frame)
##  $ Forest   : num [1:60] 1 1 1 1 1 1 1 1 1 1 ...
##  $ N        : num [1:60] 1 1 1 1 1 2 2 2 2 2 ...
##  $ P        : num [1:60] 1 1 1 1 1 1 1 1 1 1 ...
##  $ Treatment: chr [1:60] "CK" "CK" "CK" "CK" ...
##  $ Total_C  : num [1:60] 33.3 29.2 34.8 37.2 29.3 ...
##  $ Total_N  : num [1:60] 2.24 1.72 2.22 2.21 1.83 ...
##  $ TC_TN    : num [1:60] 14.9 17 15.7 16.8 16 ...
##  $ WSOC     : num [1:60] 55.6 49.9 62.6 107.3 46 ...
##  $ NH4      : num [1:60] 3.53 2.8 1.78 1.97 2.99 ...
##  $ NO3      : num [1:60] 13.3 10.2 16.6 15 14.2 ...
##  $ availP   : num [1:60] 0.309 0.275 0.71 0.615 0.553 ...
##  $ availS   : num [1:60] 48.6 41.3 32.9 29.6 33.3 ...
##  $ pH       : num [1:60] 3.7 3.66 3.56 3.58 3.64 3.66 3.48 3.55 3.68 3.5 ...
##  $ EC       : num [1:60] 119 129 151 146 154 ...
##  $ MBC      : num [1:60] 363 356 399 349 343 ...
##  $ MBN      : num [1:60] 70.4 72.7 80.9 50.9 65.8 ...
##  $ MBC_MBN  : num [1:60] 6.02 5.71 5.76 8.01 6.09 ...
##  $ ARS      : num [1:60] 49.1 46.6 52.5 53.8 52.4 ...
##  $ SARSmbc  : num [1:60] 0.135 0.131 0.132 0.154 0.153 ...
##  $ SARSt    : num [1:60] 1.47 1.59 1.51 1.45 1.79 ...
##  $ BG       : num [1:60] 244 278 312 255 326 ...
##  $ SBG      : num [1:60] 0.67 0.781 0.781 0.731 0.949 ...
##  $ NAG      : num [1:60] 216 178 372 172 239 ...
##  $ SNAG     : num [1:60] 0.595 0.501 0.932 0.492 0.696 ...
##  $ XYL      : num [1:60] 300 310 306 337 323 ...
##  $ SXYL     : num [1:60] 0.826 0.872 0.766 0.965 0.94 ...
##  $ LAP      : num [1:60] 21.6 24.8 24.6 18.8 18 ...
##  $ SLAP     : num [1:60] 0.0593 0.0695 0.0615 0.0538 0.0524 ...
##  $ ACP      : num [1:60] 18969 9279 6488 9556 10707 ...
##  $ SACP     : num [1:60] 52.2 26.1 16.2 27.4 31.2 ...
##  $ logS_C   : num [1:60] 0.618 0.602 0.616 0.624 0.611 ...
##  $ logS_N   : num [1:60] 0.712 0.723 0.662 0.759 0.713 ...
##  $ logS_P   : num [1:60] 0.395 0.42 0.451 0.435 0.427 ...
mydata$Forest<-factor(mydata$Forest,levels = c(1,2,3),labels = c("PF","RF","DF"))
mydata$Forest2<-factor(mydata$Forest,levels = c("PF","RF","DF"),labels = c("Primary Forest","Rehabilitated Forest","Disturbed Forest"))
mydata$N<-factor(mydata$N,levels = c(1,2),labels = c("N0","N150"))
mydata$P<-factor(mydata$P,levels = c(1,2),labels = c("P0","P150"))
mydata$Treatment<-factor(mydata$Treatment,levels = c("CK","N","P","NP"),labels = c("CK","N","P","NP"))

#data transform for regression
mydata.log<-data.frame(mydata[1:4],log(mydata[,c(5:12)]+1),mydata[13],
                       log(mydata[,c(14:30)]+1),mydata[31:34])

#data transform for statistical analysis
mydata.log2<-data.frame(mydata[1:4],log(mydata[,c(5:12)]),mydata[13],
                        log(mydata[,c(14:30)]),mydata[31:34])

#plot
library(vegan)
library(ggpubr)
library(readxl)
library(ggplot2)

#Fig. 1a
#PCA
library(vegan)
set.seed(123)
adonis<-adonis2(mydata[,5:17]~Forest*N*P,data = mydata,permutations = 999,method = "bray")                 
adonis
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = mydata[, 5:17] ~ Forest * N * P, data = mydata, permutations = 999, method = "bray")
##            Df SumOfSqs       R2        F Pr(>F)    
## Forest      2   1.5673  0.45191 101.6490  0.001 ***
## N           1   0.0264  0.00760   3.4199  0.035 *  
## P           1   1.7250  0.49739 223.7575  0.001 ***
## Forest:N    2   0.0258  0.00744   1.6739  0.159    
## Forest:P    2  -0.2705 -0.07800 -17.5457  1.000    
## N:P         1   0.0151  0.00436   1.9633  0.151    
## Forest:N:P  2   0.0090  0.00259   0.5834  0.651    
## Residual   48   0.3700  0.10670                    
## Total      59   3.4681  1.00000                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
data.ord<-mydata
colnames(data.ord)
##  [1] "Forest"    "N"         "P"         "Treatment" "Total_C"   "Total_N"  
##  [7] "TC_TN"     "WSOC"      "NH4"       "NO3"       "availP"    "availS"   
## [13] "pH"        "EC"        "MBC"       "MBN"       "MBC_MBN"   "ARS"      
## [19] "SARSmbc"   "SARSt"     "BG"        "SBG"       "NAG"       "SNAG"     
## [25] "XYL"       "SXYL"      "LAP"       "SLAP"      "ACP"       "SACP"     
## [31] "logS_C"    "logS_N"    "logS_P"    "Forest2"
pca<-prcomp(data.ord[,c(5:17)],scale. = TRUE)
pca
## Standard deviations (1, .., p=13):
##  [1] 2.61106181 1.57448982 1.06089326 0.92335728 0.80663187 0.56308332
##  [7] 0.54809885 0.46327968 0.32990260 0.28215176 0.21973926 0.05681361
## [13] 0.05038036
## 
## Rotation (n x k) = (13 x 13):
##                 PC1         PC2          PC3         PC4         PC5
## Total_C -0.36439584  0.03837000  0.137217479 -0.20596512  0.05867890
## Total_N -0.37210753 -0.02244927 -0.014407924  0.01181656  0.06315191
## TC_TN   -0.11432172  0.20222533  0.574874493 -0.70339382 -0.02857393
## WSOC    -0.23280945 -0.41126889 -0.112238039 -0.24623510 -0.01600838
## NH4     -0.23525969 -0.20310170 -0.264447430 -0.10875217 -0.73561753
## NO3     -0.31843461  0.23047746 -0.003701715  0.19693792  0.10081928
## availP  -0.02944763 -0.53330539  0.127331621 -0.01160790  0.46303372
## availS  -0.17094228  0.48920314 -0.018176389  0.07301420 -0.15682820
## pH       0.30569171 -0.25344614 -0.006592140 -0.22109512 -0.25774976
## EC      -0.36528912  0.04065440  0.011717626  0.17741966  0.13229570
## MBC     -0.35840990 -0.16210917 -0.018562525  0.04648275 -0.03071429
## MBN     -0.32891341 -0.23353058  0.167594344  0.18421945 -0.10376330
## MBC_MBN -0.10545405  0.15224465 -0.723037106 -0.47907047  0.32489325
##                 PC6          PC7         PC8         PC9         PC10
## Total_C  0.03606198 -0.118315994 -0.11137249 -0.35198276  0.341692996
## Total_N  0.08152112 -0.201517143 -0.11313135 -0.38670423  0.456420136
## TC_TN   -0.08308426  0.173696760  0.07011910  0.08706354 -0.155272968
## WSOC    -0.17190205 -0.556128055 -0.17551868  0.57618894  0.036580507
## NH4     -0.10139456  0.486205945 -0.13915271  0.03073294  0.103179534
## NO3     -0.38311374 -0.003469945 -0.54571942 -0.13189383 -0.530475844
## availP   0.33208727  0.479787145 -0.34354163  0.09513967  0.007172691
## availS   0.70067226 -0.067396593 -0.26346282  0.36589470  0.004477665
## pH       0.30924600 -0.301126835 -0.34837681 -0.43738973 -0.327174711
## EC      -0.09232694  0.150071956 -0.01889438  0.09023699  0.010526538
## MBC      0.19225484 -0.083195462  0.37711873 -0.13473744 -0.313920335
## MBN      0.22029948 -0.051106031  0.38604978 -0.03722634 -0.347135430
## MBC_MBN  0.08757186  0.109839186  0.14659785 -0.08204118 -0.176503570
##                PC11         PC12         PC13
## Total_C -0.02756986  0.203396841 -0.700407468
## Total_N -0.12602835 -0.272794833  0.590324834
## TC_TN    0.02579549 -0.065523555  0.203314281
## WSOC    -0.01001192  0.003265111 -0.009696621
## NH4     -0.05331246  0.005206253  0.003089812
## NO3     -0.22978057 -0.016179737  0.007576697
## availP  -0.12130498  0.037661382  0.015332711
## availS  -0.04213588  0.024395028 -0.006031454
## pH       0.34532981 -0.027943764  0.029878391
## EC       0.87957675 -0.008370397  0.049114101
## MBC     -0.09264389  0.687755329  0.239628869
## MBN     -0.09401326 -0.613832381 -0.238311408
## MBC_MBN  0.01140767 -0.164188364 -0.040713501
summary(pca)
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.6111 1.5745 1.06089 0.92336 0.80663 0.56308 0.54810
## Proportion of Variance 0.5244 0.1907 0.08658 0.06558 0.05005 0.02439 0.02311
## Cumulative Proportion  0.5244 0.7151 0.80170 0.86729 0.91734 0.94173 0.96484
##                            PC8     PC9    PC10    PC11    PC12    PC13
## Standard deviation     0.46328 0.32990 0.28215 0.21974 0.05681 0.05038
## Proportion of Variance 0.01651 0.00837 0.00612 0.00371 0.00025 0.00020
## Cumulative Proportion  0.98135 0.98972 0.99584 0.99956 0.99980 1.00000
library(ggbiplot)
ggbiplot(pca,obs.scale = 1,var.scale = 1,ellipse = TRUE) +
  scale_color_manual(name = "Treatment",values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  scale_shape_discrete(name="Forest")+
  geom_point(aes(shape=mydata$Forest,colour=mydata$Treatment),size=5)+theme_bw()+
  theme(legend.direction='vertical',legend.position='right')+
  theme(axis.text.x = element_text(size = 15,color="black"))+theme(axis.text.y = element_text(size = 15,color="black"))+
  theme(axis.title = element_text(size = 15,color="black"))+
  annotate("text",x=1.25,y=-2.5,label="Permanova analysis:\nForest: R2 = 0.45;p = 0.001\nN: R2 = 0.0076 p = 0.035\nP: R2 = 0.50; p = 0.001",color="black",size=3)

#Fig. 1b
library(ggpubr)

label = "p.signif"
method = "wilcox.test"
my_comparison<-list(c("CK","N"),c("CK","P"),c("CK","NP"),c("P","NP"))
size=17

g1<-ggboxplot(mydata,x='Treatment',y='availS',fill = 'Treatment',facet.by = 'Forest',
              add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Extractable S (mg kg-1)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g1

#ggsave("plot2/Fig.1b_availS.pdf",plot = last_plot(),units = "in",height = 5,width = 7,dpi = 600)
#Fig. 2
#ARS#
g3<-ggboxplot(mydata,x='Treatment',y='ARS',fill = 'Treatment',facet.by = 'Forest',
              add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Absolute ARS activity\n(nmol/g/h)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g3

#"BG"
g4<-ggboxplot(mydata,x='Treatment',y='BG',fill = 'Treatment',facet.by = 'Forest',
              add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Absolute BG activity\n(nmol/g/h)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g4

#"XYL"
g5<-ggboxplot(mydata,x='Treatment',y='XYL',fill = 'Treatment',facet.by = 'Forest',
              add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Absolute XYL activity\n(nmol/g/h)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g5

#"NAG"
g6<-ggboxplot(mydata.log,x='Treatment',y='NAG',fill = 'Treatment',facet.by = 'Forest',
              add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Absolute NAG activity\nLog(nmol/g/h)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g6

#"LAP"
g7<-ggboxplot(mydata,x='Treatment',y='LAP',fill = 'Treatment',facet.by = 'Forest',
              add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Absolute LAP activity\n(nmol/g/h)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g7

#"ACP"
g8<-ggboxplot(mydata.log,x='Treatment',y='ACP',fill = 'Treatment',facet.by = 'Forest',
              add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Absolute AP activity\nLog(nmol/g/h)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g8

ggarrange(g3,g4,g5,g6,g7,g8,ncol = 2,nrow = 3,widths=c(1,1),heights=c(1,1),
          common.legend = TRUE,labels = c("(a)","(b)","(c)","(d)","(e)","(f)"))

#ggsave("plot2/Fig.3_enzymeActivity2.pdf",plot = last_plot(),units = "in",height = 15,width = 14,dpi = 600)
#Fig. 3
g9<-ggboxplot(mydata.log,x='Treatment',y='SARSmbc',fill = 'Treatment',facet.by = 'Forest',
              add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Specific ARS activity\nLog(nmol ug MBC-1 h-1)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g9

g10<-ggboxplot(mydata,x='Treatment',y='logS_C',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Log(ARS)/Log(BG+XYL)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g10

#logS_N#
g11<-ggboxplot(mydata,x='Treatment',y='logS_N',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Log(ARS)/Log(NAG+LAP)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g11

#logS_P#
g12<-ggboxplot(mydata,x='Treatment',y='logS_P',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Log(ARS)/Log(AP)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g12

ggarrange(g9,g10,g11,g12,ncol = 2,nrow = 2,widths=c(1,1),heights=c(1,1),
          common.legend = TRUE,labels = c("(a)","(b)","(c)","(d)"))

#ggsave("plot2/Fig.3.pdf",plot = last_plot(),units = "in",height = 10,width = 14,dpi = 600)
####Fig. 4
#Fig. 4a
dataCON<-subset(mydata,Treatment=="CK")
dataCON.corr<-dataCON[,c(5:6,8:20)]
#controls of three forest
library(linkET)
qcorrplot(correlate(dataCON.corr,method = "spearman"), diag = FALSE,
          grid_size = 0.25) +
  geom_square() +
  geom_mark(sep = '\n',size=5,sig_level = c(0.05,0.01,0.001),sig_thres = 0.05,
            only_mark = TRUE)+
  scale_fill_gradientn(colours = RColorBrewer::brewer.pal(11, "RdYlBu"),name="rho") +
  guides(fill = guide_colorbar(title = "Spearman's rho", order = 3))

library(Hmisc)
res.cor<-rcorr(as.matrix(dataCON.corr),type = c("spearman"))
res.cor
##         Total_C Total_N  WSOC   NH4   NO3 availP availS    pH    EC   MBC   MBN
## Total_C    1.00    0.82  0.77  0.48  0.78   0.74   0.45 -0.80  0.71  0.74  0.80
## Total_N    0.82    1.00  0.62  0.40  0.91   0.57   0.75 -0.87  0.85  0.92  0.89
## WSOC       0.77    0.62  1.00  0.60  0.66   0.60   0.25 -0.68  0.61  0.61  0.61
## NH4        0.48    0.40  0.60  1.00  0.33   0.20   0.25 -0.25  0.37  0.44  0.30
## NO3        0.78    0.91  0.66  0.33  1.00   0.71   0.56 -0.87  0.82  0.81  0.85
## availP     0.74    0.57  0.60  0.20  0.71   1.00   0.05 -0.62  0.53  0.41  0.50
## availS     0.45    0.75  0.25  0.25  0.56   0.05   1.00 -0.71  0.79  0.78  0.69
## pH        -0.80   -0.87 -0.68 -0.25 -0.87  -0.62  -0.71  1.00 -0.96 -0.85 -0.87
## EC         0.71    0.85  0.61  0.37  0.82   0.53   0.79 -0.96  1.00  0.83  0.80
## MBC        0.74    0.92  0.61  0.44  0.81   0.41   0.78 -0.85  0.83  1.00  0.92
## MBN        0.80    0.89  0.61  0.30  0.85   0.50   0.69 -0.87  0.80  0.92  1.00
## MBC_MBN    0.30    0.45  0.45  0.63  0.34   0.10   0.49 -0.45  0.60  0.57  0.26
## ARS        0.70    0.62  0.66  0.79  0.63   0.51   0.25 -0.48  0.51  0.62  0.54
## SARSmbc   -0.60   -0.81 -0.53 -0.22 -0.79  -0.37  -0.80  0.90 -0.87 -0.88 -0.89
## SARSt     -0.67   -0.66 -0.64  0.01 -0.69  -0.58  -0.42  0.80 -0.65 -0.58 -0.68
##         MBC_MBN   ARS SARSmbc SARSt
## Total_C    0.30  0.70   -0.60 -0.67
## Total_N    0.45  0.62   -0.81 -0.66
## WSOC       0.45  0.66   -0.53 -0.64
## NH4        0.63  0.79   -0.22  0.01
## NO3        0.34  0.63   -0.79 -0.69
## availP     0.10  0.51   -0.37 -0.58
## availS     0.49  0.25   -0.80 -0.42
## pH        -0.45 -0.48    0.90  0.80
## EC         0.60  0.51   -0.87 -0.65
## MBC        0.57  0.62   -0.88 -0.58
## MBN        0.26  0.54   -0.89 -0.68
## MBC_MBN    1.00  0.58   -0.39 -0.09
## ARS        0.58  1.00   -0.36 -0.11
## SARSmbc   -0.39 -0.36    1.00  0.71
## SARSt     -0.09 -0.11    0.71  1.00
## 
## n= 15 
## 
## 
## P
##         Total_C Total_N WSOC   NH4    NO3    availP availS pH     EC     MBC   
## Total_C         0.0002  0.0008 0.0687 0.0006 0.0018 0.0953 0.0004 0.0028 0.0018
## Total_N 0.0002          0.0134 0.1358 0.0000 0.0261 0.0012 0.0000 0.0000 0.0000
## WSOC    0.0008  0.0134         0.0172 0.0069 0.0189 0.3760 0.0051 0.0148 0.0156
## NH4     0.0687  0.1358  0.0172        0.2265 0.4668 0.3618 0.3680 0.1728 0.1045
## NO3     0.0006  0.0000  0.0069 0.2265        0.0028 0.0297 0.0000 0.0002 0.0002
## availP  0.0018  0.0261  0.0189 0.4668 0.0028        0.8695 0.0142 0.0412 0.1283
## availS  0.0953  0.0012  0.3760 0.3618 0.0297 0.8695        0.0028 0.0005 0.0006
## pH      0.0004  0.0000  0.0051 0.3680 0.0000 0.0142 0.0028        0.0000 0.0000
## EC      0.0028  0.0000  0.0148 0.1728 0.0002 0.0412 0.0005 0.0000        0.0001
## MBC     0.0018  0.0000  0.0156 0.1045 0.0002 0.1283 0.0006 0.0000 0.0001       
## MBN     0.0003  0.0000  0.0164 0.2773 0.0000 0.0598 0.0042 0.0000 0.0004 0.0000
## MBC_MBN 0.2834  0.0895  0.0924 0.0115 0.2212 0.7134 0.0664 0.0917 0.0189 0.0272
## ARS     0.0034  0.0127  0.0073 0.0005 0.0121 0.0537 0.3688 0.0718 0.0537 0.0134
## SARSmbc 0.0181  0.0003  0.0412 0.4277 0.0005 0.1773 0.0003 0.0000 0.0000 0.0000
## SARSt   0.0065  0.0073  0.0103 0.9597 0.0045 0.0238 0.1143 0.0004 0.0087 0.0238
##         MBN    MBC_MBN ARS    SARSmbc SARSt 
## Total_C 0.0003 0.2834  0.0034 0.0181  0.0065
## Total_N 0.0000 0.0895  0.0127 0.0003  0.0073
## WSOC    0.0164 0.0924  0.0073 0.0412  0.0103
## NH4     0.2773 0.0115  0.0005 0.4277  0.9597
## NO3     0.0000 0.2212  0.0121 0.0005  0.0045
## availP  0.0598 0.7134  0.0537 0.1773  0.0238
## availS  0.0042 0.0664  0.3688 0.0003  0.1143
## pH      0.0000 0.0917  0.0718 0.0000  0.0004
## EC      0.0004 0.0189  0.0537 0.0000  0.0087
## MBC     0.0000 0.0272  0.0134 0.0000  0.0238
## MBN            0.3412  0.0380 0.0000  0.0054
## MBC_MBN 0.3412         0.0228 0.1475  0.7613
## ARS     0.0380 0.0228         0.1866  0.6851
## SARSmbc 0.0000 0.1475  0.1866         0.0028
## SARSt   0.0054 0.7613  0.6851 0.0028
library(corrplot)#corrplot()#
#corrplot(res.cor$r,type = "upper",order = "original",p.mat = res.cor$P,sig.level = 0.05, insig = "blank",tl.pos = "tp")
#corrplot(res.cor$r,add = TRUE,type = "lower",method = "number",
         #order = "original",number.cex = 0.6,diag = FALSE,tl.pos = "n",cl.pos = "n")

#Fig. 4b
dataCON<-mydata[mydata$Treatment=="CK",]

gr1<-ggplot(data = dataCON,aes(x=availS,y=ARS))+geom_point(aes(shape=Forest),size=5,color=c("#fdae61"))+
  geom_smooth(method = "lm",se=FALSE)+
  theme_bw()+ylab("Absolute ARS activity\n(nmol/g/h)")+xlab("Extractable S (mg/kg)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size),legend.position = "top")+
  annotate("text",x=42,y=25,label="R2: 0.04\n  p = 0.491",color="blue",size=5)

summary(lm(ARS~availS,data = dataCON))
## 
## Call:
## lm(formula = ARS ~ availS, data = dataCON)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.921  -7.167   1.928   4.738  13.820 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  33.9227     9.1532   3.706  0.00264 **
## availS        0.2047     0.2890   0.708  0.49140   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.854 on 13 degrees of freedom
## Multiple R-squared:  0.03714,    Adjusted R-squared:  -0.03693 
## F-statistic: 0.5014 on 1 and 13 DF,  p-value: 0.4914
gr2<-ggplot(data = dataCON,aes(x=availS,y=SARSmbc))+geom_point(aes(shape=Forest),size=5,color=c("#fdae61"))+
  geom_smooth(method = "lm",se=FALSE,color="red")+
  theme_bw()+ylab("Specific ARS activity\n(nmol/ug MBC/h)")+xlab("Extractable S (mg/kg)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size),legend.position = "top")+
  annotate("text",x=42,y=0.3,label="R2: 0.58\n  p < 0.001",color="blue",size=5)

summary(lm(SARSmbc~availS,data = dataCON))
## 
## Call:
## lm(formula = SARSmbc ~ availS, data = dataCON)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.09241 -0.03523 -0.01280  0.05062  0.09131 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.462208   0.057425   8.049 2.09e-06 ***
## availS      -0.007714   0.001813  -4.254  0.00094 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06182 on 13 degrees of freedom
## Multiple R-squared:  0.5819, Adjusted R-squared:  0.5498 
## F-statistic:  18.1 on 1 and 13 DF,  p-value: 0.0009405
library(ggpubr)
ggarrange(gr1,gr2,ncol = 1,nrow = 2,widths=c(1,1),heights=c(1,1),
          common.legend = TRUE,labels = c("(b)","(c)"))

#ggsave("plot2/Fig.4_availSenzyme.pdf",plot = last_plot(),units = "in",height = 7,width = 4,dpi = 600)

#Fig. 4c
####path analysis of dataCON####
library(lavaan)
library(semPlot)
dataCON.log<-mydata.log[mydata.log$Treatment == "CK",]

#ARS
#Fit good
model<-'
ARS~c1*MBC+c2*availS
MBC~b1*EC+b2*Total_C
availS~b3*availP+b4*EC
EC~a1*Total_C
availP~a2*Total_C

direct1:=a1*b1*c1
direct2:=a1*b4*c2
direct3:=a2*b3*c2
'
fit<-sem(model = model,data = dataCON.log,sample.nobs = 15)
summary(fit,standardized=TRUE,fit.measures=TRUE,rsq=TRUE)
## lavaan 0.6-12 ended normally after 2 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        13
## 
##   Number of observations                            15
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 8.826
##   Degrees of freedom                                 7
##   P-value (Chi-square)                           0.265
## 
## Model Test Baseline Model:
## 
##   Test statistic                               109.611
##   Degrees of freedom                                15
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.981
##   Tucker-Lewis Index (TLI)                       0.959
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)                 43.366
##   Loglikelihood unrestricted model (H1)         47.779
##                                                       
##   Akaike (AIC)                                 -60.733
##   Bayesian (BIC)                               -51.528
##   Sample-size adjusted Bayesian (BIC)          -91.216
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.132
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.361
##   P-value RMSEA <= 0.05                          0.290
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.068
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   ARS ~                                                                 
##     MBC       (c1)    0.538    0.100    5.383    0.000    0.538    1.058
##     availS    (c2)   -0.526    0.173   -3.039    0.002   -0.526   -0.597
##   MBC ~                                                                 
##     EC        (b1)    1.140    0.214    5.334    0.000    1.140    0.813
##     Total_C   (b2)    0.218    0.204    1.071    0.284    0.218    0.163
##   availS ~                                                              
##     availP    (b3)   -1.463    0.484   -3.020    0.003   -1.463   -0.552
##     EC        (b4)    0.820    0.148    5.539    0.000    0.820    1.013
##   EC ~                                                                  
##     Total_C   (a1)    0.839    0.117    7.180    0.000    0.839    0.880
##   availP ~                                                              
##     Total_C   (a2)    0.195    0.056    3.486    0.000    0.195    0.669
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .ARS               0.022    0.008    2.739    0.006    0.022    0.339
##    .MBC               0.019    0.007    2.739    0.006    0.019    0.079
##    .availS            0.027    0.010    2.739    0.006    0.027    0.328
##    .EC                0.028    0.010    2.739    0.006    0.028    0.225
##    .availP            0.007    0.002    2.739    0.006    0.007    0.552
## 
## R-Square:
##                    Estimate
##     ARS               0.661
##     MBC               0.921
##     availS            0.672
##     EC                0.775
##     availP            0.448
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     direct1           0.515    0.154    3.351    0.001    0.515    0.757
##     direct2          -0.362    0.145   -2.498    0.012   -0.362   -0.532
##     direct3           0.150    0.082    1.825    0.068    0.150    0.221
semPaths(fit,'std',layout = "tree2",sizeMan = 10,sizeLat = 10)

semPlot::semPaths(fit,'std',style="lisrel",layout = "tree2",
                  nodeLabels = c("ARS","MBC","Extractable S","EC","Extractable P","Total C"),
                  fade=FALSE,edge.label.cex = 1.75,sizeMan = 8,
                  residuals=FALSE,exoCov=FALSE)

#SARS
#Fit good
model<-'
SARSmbc~c1*MBC+c2*availS
MBC~b1*EC+b2*Total_C
availS~b3*availP+b4*EC
EC~a1*Total_C
availP~a2*Total_C

direct1:=a1*b1*c1
direct2:=a1*b4*c2
direct3:=a2*b3*c2
'
fit<-sem(model = model,data = dataCON.log,sample.nobs = 15)
summary(fit,standardized=TRUE,fit.measures=TRUE,rsq=TRUE)
## lavaan 0.6-12 ended normally after 3 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        13
## 
##   Number of observations                            15
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 8.930
##   Degrees of freedom                                 7
##   P-value (Chi-square)                           0.258
## 
## Model Test Baseline Model:
## 
##   Test statistic                               119.494
##   Degrees of freedom                                15
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.982
##   Tucker-Lewis Index (TLI)                       0.960
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)                 66.933
##   Loglikelihood unrestricted model (H1)         71.398
##                                                       
##   Akaike (AIC)                                -107.866
##   Bayesian (BIC)                               -98.661
##   Sample-size adjusted Bayesian (BIC)         -138.349
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.136
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.363
##   P-value RMSEA <= 0.05                          0.282
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.037
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   SARSmbc ~                                                             
##     MBC       (c1)   -0.073    0.021   -3.516    0.000   -0.073   -0.512
##     availS    (c2)   -0.119    0.036   -3.323    0.001   -0.119   -0.484
##   MBC ~                                                                 
##     EC        (b1)    1.140    0.214    5.334    0.000    1.140    0.813
##     Total_C   (b2)    0.218    0.204    1.071    0.284    0.218    0.163
##   availS ~                                                              
##     availP    (b3)   -1.463    0.484   -3.020    0.003   -1.463   -0.552
##     EC        (b4)    0.820    0.148    5.539    0.000    0.820    1.013
##   EC ~                                                                  
##     Total_C   (a1)    0.839    0.117    7.180    0.000    0.839    0.880
##   availP ~                                                              
##     Total_C   (a2)    0.195    0.056    3.486    0.000    0.195    0.669
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .SARSmbc           0.001    0.000    2.739    0.006    0.001    0.186
##    .MBC               0.019    0.007    2.739    0.006    0.019    0.079
##    .availS            0.027    0.010    2.739    0.006    0.027    0.328
##    .EC                0.028    0.010    2.739    0.006    0.028    0.225
##    .availP            0.007    0.002    2.739    0.006    0.007    0.552
## 
## R-Square:
##                    Estimate
##     SARSmbc           0.814
##     MBC               0.921
##     availS            0.672
##     EC                0.775
##     availP            0.448
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     direct1          -0.070    0.026   -2.717    0.007   -0.070   -0.366
##     direct2          -0.082    0.031   -2.649    0.008   -0.082   -0.431
##     direct3           0.034    0.018    1.881    0.060    0.034    0.179
semPaths(fit,'std',layout = "tree2",sizeMan = 10,sizeLat = 10)

semPlot::semPaths(fit,'std',style="lisrel",layout = "tree2",
                  nodeLabels = c("SARS","MBC","Extractable S","EC","Extractable P","Total C"),
                  fade=FALSE,edge.label.cex = 1.75,sizeMan = 8,
                  residuals=FALSE,exoCov=FALSE)

##Fig. 5
###Primary forest
library(psych)
library(tidyverse)
dataMF<-mydata[mydata$Forest == "PF",]
colnames(dataMF)
##  [1] "Forest"    "N"         "P"         "Treatment" "Total_C"   "Total_N"  
##  [7] "TC_TN"     "WSOC"      "NH4"       "NO3"       "availP"    "availS"   
## [13] "pH"        "EC"        "MBC"       "MBN"       "MBC_MBN"   "ARS"      
## [19] "SARSmbc"   "SARSt"     "BG"        "SBG"       "NAG"       "SNAG"     
## [25] "XYL"       "SXYL"      "LAP"       "SLAP"      "ACP"       "SACP"     
## [31] "logS_C"    "logS_N"    "logS_P"    "Forest2"
corr.tmp1<-dataMF[,c("Total_C","Total_N","TC_TN","WSOC","NH4","NO3","availP","availS",
                     "pH","EC","MBC","MBN","MBC_MBN")]
corr.tmp2<-dataMF[,c("ARS","SARSmbc","SARSt","logS_C","logS_N","logS_P")]

spman.d12 = corr.test(corr.tmp1, corr.tmp2,use="pairwise",method="spearman",adjust="fdr",alpha=.05,ci=FALSE)

r<-data.frame(spman.d12$r)
r[is.na(r)] <-0
r$X<-row.names(r)

r.long <- gather(r, Y, r, 1 : ncol(r)-1,  factor_key=TRUE)
p<-data.frame(spman.d12$p)
p$X<-row.names(p)
p.long <- gather(p, Y, p, 1 : ncol(p)-1, factor_key=TRUE)
cor.out<-cbind(r.long,p.long$p)
cor.out$r <- round(as.numeric(cor.out$r), 2)
str(cor.out$r)
##  num [1:78] -0.12 0.01 -0.09 0.39 0.07 -0.52 0.69 -0.69 0.26 0.13 ...
cor.out$X<-factor(cor.out$X,levels = c("Total_C","Total_N","TC_TN","WSOC","NH4","NO3","availP","availS",
                                       "pH","EC","MBC","MBN","MBC_MBN"))
cor.out$X<-factor(cor.out$X,labels = c("Total C","Total N","Total C/total N","WSOC","NH4+","NO3-","extractable P","extractable S",
                                       "pH","EC","MBC","MBN","MBC/MBN"))

cor.out$Y<-factor(cor.out$Y,levels = c("ARS","SARSmbc","SARSt","logS_C","logS_N","logS_P"))
cor.out$Y<-factor(cor.out$Y,labels = c("ARS","SARSmbc","SARStoc","LogS_C","LogS_N","LogS_P"))

ggplot(cor.out, aes(Y,X)) + 
  geom_tile(aes(fill = r), size=1)+
  geom_text(aes(label = r), size = 1.5)+
  scale_fill_gradient(guide = "legend", high='Salmon', low='CornflowerBlue')+
  theme(axis.title= element_blank(),
        axis.text.x=element_text(colour="black", size=8, face="bold",angle = 20),
        axis.text.y=element_text(colour="black", size=8, face="bold"))

cor.out.sig <- cor.out
cor.out.sig$r[cor.out.sig$p > 0.05] <- 0
cor.out$Forest<-"PF"
cor.out.sig$Forest<-"PF"

ggplot(cor.out.sig, aes(Y,X)) + 
  geom_tile(aes(fill = r), size=1)+
  geom_text(data= cor.out.sig[cor.out.sig$r!=0,], aes(label = r), size = 5, font="bold")+
  scale_fill_gradient2(guide = "legend", high='#2c7bb6',mid="white", low='#d7191c')+
  theme_bw()+
  theme(axis.title= element_blank(),
        axis.text.x=element_text(colour="black", size=15, face="bold",angle = 30, hjust=1),
        axis.text.y=element_text(colour="black", size=15, face="bold"))

cor.out.PF<-cor.out
cor.out.sig.PF<-cor.out.sig


###Rehabilitated forest
library(psych)
library(tidyverse)
dataRF<-mydata[mydata$Forest=="RF",]
colnames(dataRF)
##  [1] "Forest"    "N"         "P"         "Treatment" "Total_C"   "Total_N"  
##  [7] "TC_TN"     "WSOC"      "NH4"       "NO3"       "availP"    "availS"   
## [13] "pH"        "EC"        "MBC"       "MBN"       "MBC_MBN"   "ARS"      
## [19] "SARSmbc"   "SARSt"     "BG"        "SBG"       "NAG"       "SNAG"     
## [25] "XYL"       "SXYL"      "LAP"       "SLAP"      "ACP"       "SACP"     
## [31] "logS_C"    "logS_N"    "logS_P"    "Forest2"
corr.tmp1<-dataRF[,c("Total_C","Total_N","TC_TN","WSOC","NH4","NO3","availP","availS",
                     "pH","EC","MBC","MBN","MBC_MBN")]
corr.tmp2<-dataRF[,c("ARS","SARSmbc","SARSt","logS_C","logS_N","logS_P")]

spman.d12 = corr.test(corr.tmp1, corr.tmp2,use="pairwise",method="spearman",adjust="fdr",alpha=.05,ci=FALSE)

r<-data.frame(spman.d12$r)
r[is.na(r)] <-0
r$X<-row.names(r)

r.long <- gather(r, Y, r, 1 : ncol(r)-1,  factor_key=TRUE)
p<-data.frame(spman.d12$p)
p$X<-row.names(p)
p.long <- gather(p, Y, p, 1 : ncol(p)-1, factor_key=TRUE)
cor.out<-cbind(r.long,p.long$p)
cor.out$r <- round(as.numeric(cor.out$r), 2)
str(cor.out$r)
##  num [1:78] 0.15 0.3 0.21 -0.18 -0.01 -0.01 0.08 -0.24 0.13 -0.15 ...
cor.out$X<-factor(cor.out$X,levels = c("Total_C","Total_N","TC_TN","WSOC","NH4","NO3","availP","availS",
                                       "pH","EC","MBC","MBN","MBC_MBN"))
cor.out$X<-factor(cor.out$X,labels = c("Total C","Total N","Total C/total N","WSOC","NH4+","NO3-","extractable P","extractable S",
                                       "pH","EC","MBC","MBN","MBC/MBN"))

cor.out$Y<-factor(cor.out$Y,levels = c("ARS","SARSmbc","SARSt","logS_C","logS_N","logS_P"))
cor.out$Y<-factor(cor.out$Y,labels = c("ARS","SARSmbc","SARStoc","LogS_C","LogS_N","LogS_P"))

ggplot(cor.out, aes(Y,X)) + 
  geom_tile(aes(fill = r), size=1)+
  geom_text(aes(label = r), size = 1.5)+
  scale_fill_gradient(guide = "legend", high='Salmon', low='CornflowerBlue')+
  theme(axis.title= element_blank(),
        axis.text.x=element_text(colour="black", size=8, face="bold",angle = 20),
        axis.text.y=element_text(colour="black", size=8, face="bold"))

cor.out.sig <- cor.out
cor.out.sig$r[cor.out.sig$p > 0.05] <- 0
cor.out$Forest<-"RF"
cor.out.sig$Forest<-"RF"

ggplot(cor.out.sig, aes(Y,X)) + 
  geom_tile(aes(fill = r), size=1)+
  geom_text(data= cor.out.sig[cor.out.sig$r!=0,], aes(label = r), size = 5, font="bold")+
  scale_fill_gradient2(guide = "legend", high='#2c7bb6',mid="white", low='#d7191c')+
  theme_bw()+
  theme(axis.title= element_blank(),
        axis.text.x=element_text(colour="black", size=15, face="bold",angle = 30, hjust=1),
        axis.text.y=element_text(colour="black", size=15, face="bold"))

cor.out.RF<-cor.out
cor.out.sig.RF<-cor.out.sig


###Disturbed forest
library(psych)
library(tidyverse)
dataDF<-mydata[mydata$Forest == "DF",]
colnames(dataDF)
##  [1] "Forest"    "N"         "P"         "Treatment" "Total_C"   "Total_N"  
##  [7] "TC_TN"     "WSOC"      "NH4"       "NO3"       "availP"    "availS"   
## [13] "pH"        "EC"        "MBC"       "MBN"       "MBC_MBN"   "ARS"      
## [19] "SARSmbc"   "SARSt"     "BG"        "SBG"       "NAG"       "SNAG"     
## [25] "XYL"       "SXYL"      "LAP"       "SLAP"      "ACP"       "SACP"     
## [31] "logS_C"    "logS_N"    "logS_P"    "Forest2"
corr.tmp1<-dataDF[,c("Total_C","Total_N","TC_TN","WSOC","NH4","NO3","availP","availS",
                     "pH","EC","MBC","MBN","MBC_MBN")]
corr.tmp2<-dataDF[,c("ARS","SARSmbc","SARSt","logS_C","logS_N","logS_P")]

spman.d12 = corr.test(corr.tmp1, corr.tmp2,use="pairwise",method="spearman",adjust="fdr",alpha=.05,ci=FALSE)

r<-data.frame(spman.d12$r)
r[is.na(r)] <-0
r$X<-row.names(r)

r.long <- gather(r, Y, r, 1 : ncol(r)-1,  factor_key=TRUE)
p<-data.frame(spman.d12$p)
p$X<-row.names(p)
p.long <- gather(p, Y, p, 1 : ncol(p)-1, factor_key=TRUE)
cor.out<-cbind(r.long,p.long$p)
cor.out$r <- round(as.numeric(cor.out$r), 2)
str(cor.out$r)
##  num [1:78] 0.03 0.23 0.1 0.58 0.23 -0.41 0.72 -0.68 0.4 -0.25 ...
cor.out$X<-factor(cor.out$X,levels = c("Total_C","Total_N","TC_TN","WSOC","NH4","NO3","availP","availS",
                                       "pH","EC","MBC","MBN","MBC_MBN"))
cor.out$X<-factor(cor.out$X,labels = c("Total C","Total N","Total C/total N","WSOC","NH4+","NO3-","extractable P","extractable S",
                                       "pH","EC","MBC","MBN","MBC/MBN"))

cor.out$Y<-factor(cor.out$Y,levels = c("ARS","SARSmbc","SARSt","logS_C","logS_N","logS_P"))
cor.out$Y<-factor(cor.out$Y,labels = c("ARS","SARSmbc","SARStoc","LogS_C","LogS_N","LogS_P"))

ggplot(cor.out, aes(Y,X)) + 
  geom_tile(aes(fill = r), size=1)+
  geom_text(aes(label = r), size = 1.5)+
  scale_fill_gradient(guide = "legend", high='Salmon', low='CornflowerBlue')+
  theme(axis.title= element_blank(),
        axis.text.x=element_text(colour="black", size=8, face="bold",angle = 20),
        axis.text.y=element_text(colour="black", size=8, face="bold"))

cor.out.sig <- cor.out
cor.out.sig$r[cor.out.sig$p > 0.05] <- 0
cor.out$Forest<-"DF"
cor.out.sig$Forest<-"DF"

ggplot(cor.out.sig, aes(Y,X)) + 
  geom_tile(aes(fill = r), size=1)+
  geom_text(data= cor.out.sig[cor.out.sig$r!=0,], aes(label = r), size = 5, font="bold")+
  scale_fill_gradient2(guide = "legend", high='#2c7bb6',mid="white", low='#d7191c')+
  theme_bw()+
  theme(axis.title= element_blank(),
        axis.text.x=element_text(colour="black", size=15, face="bold",angle = 30, hjust=1),
        axis.text.y=element_text(colour="black", size=15, face="bold"))

cor.out.DF<-cor.out
cor.out.sig.DF<-cor.out.sig


colnames(cor.out.PF) == colnames(cor.out.RF)
## [1] TRUE TRUE TRUE TRUE TRUE
colnames(cor.out.RF) == colnames(cor.out.DF)
## [1] TRUE TRUE TRUE TRUE TRUE
cor.out.all<-rbind(cor.out.PF,cor.out.RF,cor.out.DF)
cor.out.all$Forest<-factor(cor.out.all$Forest,levels = c("PF","RF","DF"),
                           labels = c("PF-Primary Forest","RF-Rehablitated Forest","DF-Disturbed Forest"))

colnames(cor.out.sig.PF) == colnames(cor.out.sig.RF)
## [1] TRUE TRUE TRUE TRUE TRUE
colnames(cor.out.sig.RF) == colnames(cor.out.sig.DF)
## [1] TRUE TRUE TRUE TRUE TRUE
cor.out.sig.all<-rbind(cor.out.sig.PF,cor.out.sig.RF,cor.out.sig.DF)
cor.out.sig.all$Forest<-factor(cor.out.sig.all$Forest,levels = c("PF","RF","DF"),
                               labels = c("PF-Primary Forest","RF-Rehablitated Forest","DF-Disturbed Forest"))

ggplot(cor.out.sig.all, aes(Y,X)) + 
  geom_tile(aes(fill = r), size=1)+
  geom_text(data= cor.out.sig.all[cor.out.sig.all$r!=0,], aes(label = r), size = 5, font="bold")+
  scale_fill_gradient2(guide = "legend", high='#2c7bb6',mid="white", low='#d7191c')+
  theme_bw()+facet_grid(~Forest)+guides(fill = guide_legend(title = 'rho'))+
  theme(axis.title= element_blank(),
        axis.text.x=element_text(colour="black", size=15, face="bold",angle = 30, hjust=1),
        axis.text.y=element_text(colour="black", size=15, face="bold"),
        strip.text.x = element_text(size = 15,colour = "black"),
        legend.key.size = unit(25,"pt"))

#ggsave("plot2/Fig.5_spearmanCorr.pdf",plot = last_plot(),units = "in",height = 6,width = 12,dpi = 600)

#Fig. 5b-c
####Data MF####
dataMF.log<-subset(mydata.log,Forest=="PF")

####stepwise regression analysis####
#ARS~
data.step<- dataMF.log[,5:18]
library(olsrr)
model<-ARS~.
fit<-lm(formula=model,data = data.step)
ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method    
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. Total_C 
## 2. Total_N 
## 3. TC_TN 
## 4. WSOC 
## 5. NH4 
## 6. NO3 
## 7. availP 
## 8. availS 
## 9. pH 
## 10. EC 
## 11. MBC 
## 12. MBN 
## 13. MBC_MBN 
## 
## We are selecting variables based on p value...
## 
## 
## Forward Selection: Step 1 
## 
## - availP 
## 
##                         Model Summary                         
## -------------------------------------------------------------
## R                       0.820       RMSE               0.156 
## R-Squared               0.672       Coef. Var          3.738 
## Adj. R-Squared          0.654       MSE                0.024 
## Pred R-Squared          0.592       MAE                0.110 
## -------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression      0.899         1          0.899    36.948    0.0000 
## Residual        0.438        18          0.024                     
## Total           1.338        19                                    
## -------------------------------------------------------------------
## 
##                                  Parameter Estimates                                  
## -------------------------------------------------------------------------------------
##       model     Beta    Std. Error    Std. Beta      t        Sig     lower    upper 
## -------------------------------------------------------------------------------------
## (Intercept)    3.904         0.057                 69.041    0.000    3.785    4.023 
##      availP    0.085         0.014        0.820     6.078    0.000    0.056    0.115 
## -------------------------------------------------------------------------------------
## 
## 
## 
## No more variables to be added.
## 
## Variables Entered: 
## 
## + availP 
## 
## 
## Final Model Output 
## ------------------
## 
##                         Model Summary                         
## -------------------------------------------------------------
## R                       0.820       RMSE               0.156 
## R-Squared               0.672       Coef. Var          3.738 
## Adj. R-Squared          0.654       MSE                0.024 
## Pred R-Squared          0.592       MAE                0.110 
## -------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression      0.899         1          0.899    36.948    0.0000 
## Residual        0.438        18          0.024                     
## Total           1.338        19                                    
## -------------------------------------------------------------------
## 
##                                  Parameter Estimates                                  
## -------------------------------------------------------------------------------------
##       model     Beta    Std. Error    Std. Beta      t        Sig     lower    upper 
## -------------------------------------------------------------------------------------
## (Intercept)    3.904         0.057                 69.041    0.000    3.785    4.023 
##      availP    0.085         0.014        0.820     6.078    0.000    0.056    0.115 
## -------------------------------------------------------------------------------------
## 
##                             Selection Summary                             
## -------------------------------------------------------------------------
##         Variable                  Adj.                                       
## Step    Entered     R-Square    R-Square     C(p)        AIC        RMSE     
## -------------------------------------------------------------------------
##    1    availP        0.6724      0.6542    39.1403    -13.6600    0.1560    
## -------------------------------------------------------------------------
test<-ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method    
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. Total_C 
## 2. Total_N 
## 3. TC_TN 
## 4. WSOC 
## 5. NH4 
## 6. NO3 
## 7. availP 
## 8. availS 
## 9. pH 
## 10. EC 
## 11. MBC 
## 12. MBN 
## 13. MBC_MBN 
## 
## We are selecting variables based on p value...
## 
## 
## Forward Selection: Step 1 
## 
## - availP 
## 
##                         Model Summary                         
## -------------------------------------------------------------
## R                       0.820       RMSE               0.156 
## R-Squared               0.672       Coef. Var          3.738 
## Adj. R-Squared          0.654       MSE                0.024 
## Pred R-Squared          0.592       MAE                0.110 
## -------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression      0.899         1          0.899    36.948    0.0000 
## Residual        0.438        18          0.024                     
## Total           1.338        19                                    
## -------------------------------------------------------------------
## 
##                                  Parameter Estimates                                  
## -------------------------------------------------------------------------------------
##       model     Beta    Std. Error    Std. Beta      t        Sig     lower    upper 
## -------------------------------------------------------------------------------------
## (Intercept)    3.904         0.057                 69.041    0.000    3.785    4.023 
##      availP    0.085         0.014        0.820     6.078    0.000    0.056    0.115 
## -------------------------------------------------------------------------------------
## 
## 
## 
## No more variables to be added.
## 
## Variables Entered: 
## 
## + availP 
## 
## 
## Final Model Output 
## ------------------
## 
##                         Model Summary                         
## -------------------------------------------------------------
## R                       0.820       RMSE               0.156 
## R-Squared               0.672       Coef. Var          3.738 
## Adj. R-Squared          0.654       MSE                0.024 
## Pred R-Squared          0.592       MAE                0.110 
## -------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression      0.899         1          0.899    36.948    0.0000 
## Residual        0.438        18          0.024                     
## Total           1.338        19                                    
## -------------------------------------------------------------------
## 
##                                  Parameter Estimates                                  
## -------------------------------------------------------------------------------------
##       model     Beta    Std. Error    Std. Beta      t        Sig     lower    upper 
## -------------------------------------------------------------------------------------
## (Intercept)    3.904         0.057                 69.041    0.000    3.785    4.023 
##      availP    0.085         0.014        0.820     6.078    0.000    0.056    0.115 
## -------------------------------------------------------------------------------------
result.stepre.ars.MF<-data.frame(test$predictors,test$adjr)
colnames(result.stepre.ars.MF)<-c("prd","adjr")
#View(result.stepre.ars.MF)
result.stepre.ars.MF$prd2<-c("extractP")
result.stepre.ars.MF$adjp<-c("<0.001")
result.stepre.ars.MF$adjsig<-c("***")
result.stepre.ars.MF$adjr2<-c(0.6542)
library(ggplot2)
p1<-ggplot()+geom_col(data=result.stepre.ars.MF,aes(x=prd2,y=adjr),fill="#fdae61",color=NA,width = 0.5)+
  labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
  theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
                   axis.title = element_text(size = 15,colour = "black"))+
  geom_text(data = result.stepre.ars.MF,aes(x=prd2,y=adjr,label=adjsig),nudge_y = 0.05,size=7)
p1

#SARSm~
set.seed(123)
data.step<- dataMF.log[,c(5:17,19)]
library(olsrr)
model<-SARSmbc~.
fit<-lm(formula=model,data = data.step)
ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method    
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. Total_C 
## 2. Total_N 
## 3. TC_TN 
## 4. WSOC 
## 5. NH4 
## 6. NO3 
## 7. availP 
## 8. availS 
## 9. pH 
## 10. EC 
## 11. MBC 
## 12. MBN 
## 13. MBC_MBN 
## 
## We are selecting variables based on p value...
## 
## 
## Forward Selection: Step 1 
## 
## - MBC_MBN 
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.557       RMSE                0.031 
## R-Squared               0.310       Coef. Var          18.784 
## Adj. R-Squared          0.272       MSE                 0.001 
## Pred R-Squared          0.001       MAE                 0.024 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.008         1          0.008    8.094    0.0107 
## Residual        0.017        18          0.001                    
## Total           0.025        19                                   
## ------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     0.469         0.107                  4.391    0.000     0.245     0.694 
##     MBC_MBN    -0.166         0.058       -0.557    -2.845    0.011    -0.288    -0.043 
## ----------------------------------------------------------------------------------------
## 
## 
## 
## Forward Selection: Step 2 
## 
## - MBC 
## 
##                          Model Summary                          
## ---------------------------------------------------------------
## R                        0.671       RMSE                0.029 
## R-Squared                0.450       Coef. Var          17.260 
## Adj. R-Squared           0.385       MSE                 0.001 
## Pred R-Squared          -0.006       MAE                 0.023 
## ---------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.011         2          0.006    6.952    0.0062 
## Residual        0.014        17          0.001                    
## Total           0.025        19                                   
## ------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     0.904         0.231                  3.910    0.001     0.416     1.392 
##     MBC_MBN    -0.183         0.054       -0.614    -3.373    0.004    -0.297    -0.068 
##         MBC    -0.069         0.033       -0.378    -2.078    0.053    -0.138     0.001 
## ----------------------------------------------------------------------------------------
## 
## 
## 
## Forward Selection: Step 3 
## 
## - availP 
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.769       RMSE                0.025 
## R-Squared               0.591       Coef. Var          15.345 
## Adj. R-Squared          0.514       MSE                 0.001 
## Pred R-Squared          0.150       MAE                 0.018 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.015         3          0.005      7.7    0.0021 
## Residual        0.010        16          0.001                    
## Total           0.025        19                                   
## ------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     0.944         0.206                  4.575    0.000     0.507     1.381 
##     MBC_MBN    -0.043         0.077       -0.144    -0.559    0.584    -0.205     0.120 
##         MBC    -0.124         0.038       -0.684    -3.293    0.005    -0.204    -0.044 
##      availP     0.010         0.004        0.706     2.347    0.032     0.001     0.019 
## ----------------------------------------------------------------------------------------
## 
## 
## 
## No more variables to be added.
## 
## Variables Entered: 
## 
## + MBC_MBN 
## + MBC 
## + availP 
## 
## 
## Final Model Output 
## ------------------
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.769       RMSE                0.025 
## R-Squared               0.591       Coef. Var          15.345 
## Adj. R-Squared          0.514       MSE                 0.001 
## Pred R-Squared          0.150       MAE                 0.018 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.015         3          0.005      7.7    0.0021 
## Residual        0.010        16          0.001                    
## Total           0.025        19                                   
## ------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     0.944         0.206                  4.575    0.000     0.507     1.381 
##     MBC_MBN    -0.043         0.077       -0.144    -0.559    0.584    -0.205     0.120 
##         MBC    -0.124         0.038       -0.684    -3.293    0.005    -0.204    -0.044 
##      availP     0.010         0.004        0.706     2.347    0.032     0.001     0.019 
## ----------------------------------------------------------------------------------------
## 
##                             Selection Summary                             
## -------------------------------------------------------------------------
##         Variable                  Adj.                                       
## Step    Entered     R-Square    R-Square     C(p)        AIC        RMSE     
## -------------------------------------------------------------------------
##    1    MBC_MBN       0.3102      0.2719    71.0093    -78.1093    0.0311    
##    2    MBC           0.4499      0.3852    55.3836    -80.6366    0.0286    
##    3    availP        0.5908      0.5141    39.6168    -84.5526    0.0254    
## -------------------------------------------------------------------------
test<-ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method    
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. Total_C 
## 2. Total_N 
## 3. TC_TN 
## 4. WSOC 
## 5. NH4 
## 6. NO3 
## 7. availP 
## 8. availS 
## 9. pH 
## 10. EC 
## 11. MBC 
## 12. MBN 
## 13. MBC_MBN 
## 
## We are selecting variables based on p value...
## 
## 
## Forward Selection: Step 1 
## 
## - MBC_MBN 
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.557       RMSE                0.031 
## R-Squared               0.310       Coef. Var          18.784 
## Adj. R-Squared          0.272       MSE                 0.001 
## Pred R-Squared          0.001       MAE                 0.024 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.008         1          0.008    8.094    0.0107 
## Residual        0.017        18          0.001                    
## Total           0.025        19                                   
## ------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     0.469         0.107                  4.391    0.000     0.245     0.694 
##     MBC_MBN    -0.166         0.058       -0.557    -2.845    0.011    -0.288    -0.043 
## ----------------------------------------------------------------------------------------
## 
## 
## 
## Forward Selection: Step 2 
## 
## - MBC 
## 
##                          Model Summary                          
## ---------------------------------------------------------------
## R                        0.671       RMSE                0.029 
## R-Squared                0.450       Coef. Var          17.260 
## Adj. R-Squared           0.385       MSE                 0.001 
## Pred R-Squared          -0.006       MAE                 0.023 
## ---------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.011         2          0.006    6.952    0.0062 
## Residual        0.014        17          0.001                    
## Total           0.025        19                                   
## ------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     0.904         0.231                  3.910    0.001     0.416     1.392 
##     MBC_MBN    -0.183         0.054       -0.614    -3.373    0.004    -0.297    -0.068 
##         MBC    -0.069         0.033       -0.378    -2.078    0.053    -0.138     0.001 
## ----------------------------------------------------------------------------------------
## 
## 
## 
## Forward Selection: Step 3 
## 
## - availP 
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.769       RMSE                0.025 
## R-Squared               0.591       Coef. Var          15.345 
## Adj. R-Squared          0.514       MSE                 0.001 
## Pred R-Squared          0.150       MAE                 0.018 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.015         3          0.005      7.7    0.0021 
## Residual        0.010        16          0.001                    
## Total           0.025        19                                   
## ------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     0.944         0.206                  4.575    0.000     0.507     1.381 
##     MBC_MBN    -0.043         0.077       -0.144    -0.559    0.584    -0.205     0.120 
##         MBC    -0.124         0.038       -0.684    -3.293    0.005    -0.204    -0.044 
##      availP     0.010         0.004        0.706     2.347    0.032     0.001     0.019 
## ----------------------------------------------------------------------------------------
## 
## 
## 
## No more variables to be added.
## 
## Variables Entered: 
## 
## + MBC_MBN 
## + MBC 
## + availP 
## 
## 
## Final Model Output 
## ------------------
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.769       RMSE                0.025 
## R-Squared               0.591       Coef. Var          15.345 
## Adj. R-Squared          0.514       MSE                 0.001 
## Pred R-Squared          0.150       MAE                 0.018 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.015         3          0.005      7.7    0.0021 
## Residual        0.010        16          0.001                    
## Total           0.025        19                                   
## ------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     0.944         0.206                  4.575    0.000     0.507     1.381 
##     MBC_MBN    -0.043         0.077       -0.144    -0.559    0.584    -0.205     0.120 
##         MBC    -0.124         0.038       -0.684    -3.293    0.005    -0.204    -0.044 
##      availP     0.010         0.004        0.706     2.347    0.032     0.001     0.019 
## ----------------------------------------------------------------------------------------
plot(test)

result.stepre.sars.MF<-data.frame(test$predictors,test$adjr)
colnames(result.stepre.sars.MF)<-c("prd","adjr")
#View(result.stepre.sars.MF)
result.stepre.sars.MF$prd2<-c("MBC_MBN","MBC","extractP")
result.stepre.sars.MF$adjp<-c("0.584","0.005","0.032")
result.stepre.sars.MF$adjsig<-c(".","***","*")
result.stepre.sars.MF$adjr2<-c(0.2719,0.1133,0.1289)
library(ggplot2)
p2<-ggplot()+geom_col(data=result.stepre.sars.MF,aes(x=prd2,y=adjr2),fill="#fdae61",color=NA,width = 0.5)+
  labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
  theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
                   axis.title = element_text(size = 15,colour = "black"))+
  geom_text(data = result.stepre.sars.MF,aes(x=prd2,y=adjr2,label=adjsig),nudge_y = 0.05,size=7)
p2

library(ggpubr)
ggarrange(p1,p2,ncol = 1,nrow = 2,widths=c(1,1),heights=c(1,1),labels = c("(b)","(c)"))

#ggsave("plot2/Fig.2_stepwiseControls.pdf",plot = last_plot(),units = "in",height = 12,width = 6,dpi = 600)

result.stepre.ars.MF$enzyme<-"ARS-PF"
result.stepre.sars.MF$enzyme<-"SARS-PF"
result.stepre.MF<-rbind(result.stepre.ars.MF,result.stepre.sars.MF)
result.stepre.MF$forest<-"PF"
#View(result.stepre.MF)

library(ggplot2)
p3<-ggplot()+geom_col(data=result.stepre.MF,aes(x=prd2,y=adjr2),fill="#fdae61",color=NA,width = 0.5)+
  labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
  theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
                   axis.title = element_text(size = 15,colour = "black"),
                   axis.text.x = element_text(angle = 15,hjust = 1),
                   strip.text.x = element_text(size = 15,colour = "black"))+
  geom_text(data = result.stepre.MF,aes(x=prd2,y=adjr,label=adjsig),nudge_y = 0.05,size=7)+
  facet_grid(~enzyme,scales = "free")
p3

####Data RF####
dataRF.log<-subset(mydata.log,Forest=="RF")

####stepwise regression analysis####
#ARS~
data.step<- dataRF.log[,5:18]
library(olsrr)
model<-ARS~.
fit<-lm(formula=model,data = data.step)
ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method    
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. Total_C 
## 2. Total_N 
## 3. TC_TN 
## 4. WSOC 
## 5. NH4 
## 6. NO3 
## 7. availP 
## 8. availS 
## 9. pH 
## 10. EC 
## 11. MBC 
## 12. MBN 
## 13. MBC_MBN 
## 
## We are selecting variables based on p value...
## 
## 
## Forward Selection: Step 1 
## 
## - MBC 
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                        0.378       RMSE               0.168 
## R-Squared                0.143       Coef. Var          4.617 
## Adj. R-Squared           0.096       MSE                0.028 
## Pred R-Squared          -0.024       MAE                0.139 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.085         1          0.085    3.009    0.0999 
## Residual        0.507        18          0.028                    
## Total           0.592        19                                   
## ------------------------------------------------------------------
## 
##                                  Parameter Estimates                                  
## -------------------------------------------------------------------------------------
##       model     Beta    Std. Error    Std. Beta      t       Sig      lower    upper 
## -------------------------------------------------------------------------------------
## (Intercept)    2.349         0.743                 3.162    0.005     0.788    3.910 
##         MBC    0.264         0.152        0.378    1.735    0.100    -0.056    0.584 
## -------------------------------------------------------------------------------------
## 
## 
## 
## No more variables to be added.
## 
## Variables Entered: 
## 
## + MBC 
## 
## 
## Final Model Output 
## ------------------
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                        0.378       RMSE               0.168 
## R-Squared                0.143       Coef. Var          4.617 
## Adj. R-Squared           0.096       MSE                0.028 
## Pred R-Squared          -0.024       MAE                0.139 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.085         1          0.085    3.009    0.0999 
## Residual        0.507        18          0.028                    
## Total           0.592        19                                   
## ------------------------------------------------------------------
## 
##                                  Parameter Estimates                                  
## -------------------------------------------------------------------------------------
##       model     Beta    Std. Error    Std. Beta      t       Sig      lower    upper 
## -------------------------------------------------------------------------------------
## (Intercept)    2.349         0.743                 3.162    0.005     0.788    3.910 
##         MBC    0.264         0.152        0.378    1.735    0.100    -0.056    0.584 
## -------------------------------------------------------------------------------------
## 
##                            Selection Summary                             
## ------------------------------------------------------------------------
##         Variable                  Adj.                                      
## Step    Entered     R-Square    R-Square     C(p)       AIC        RMSE     
## ------------------------------------------------------------------------
##    1    MBC           0.1432      0.0956    9.1247    -10.7299    0.1679    
## ------------------------------------------------------------------------
test<-ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method    
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. Total_C 
## 2. Total_N 
## 3. TC_TN 
## 4. WSOC 
## 5. NH4 
## 6. NO3 
## 7. availP 
## 8. availS 
## 9. pH 
## 10. EC 
## 11. MBC 
## 12. MBN 
## 13. MBC_MBN 
## 
## We are selecting variables based on p value...
## 
## 
## Forward Selection: Step 1 
## 
## - MBC 
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                        0.378       RMSE               0.168 
## R-Squared                0.143       Coef. Var          4.617 
## Adj. R-Squared           0.096       MSE                0.028 
## Pred R-Squared          -0.024       MAE                0.139 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.085         1          0.085    3.009    0.0999 
## Residual        0.507        18          0.028                    
## Total           0.592        19                                   
## ------------------------------------------------------------------
## 
##                                  Parameter Estimates                                  
## -------------------------------------------------------------------------------------
##       model     Beta    Std. Error    Std. Beta      t       Sig      lower    upper 
## -------------------------------------------------------------------------------------
## (Intercept)    2.349         0.743                 3.162    0.005     0.788    3.910 
##         MBC    0.264         0.152        0.378    1.735    0.100    -0.056    0.584 
## -------------------------------------------------------------------------------------
## 
## 
## 
## No more variables to be added.
## 
## Variables Entered: 
## 
## + MBC 
## 
## 
## Final Model Output 
## ------------------
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                        0.378       RMSE               0.168 
## R-Squared                0.143       Coef. Var          4.617 
## Adj. R-Squared           0.096       MSE                0.028 
## Pred R-Squared          -0.024       MAE                0.139 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.085         1          0.085    3.009    0.0999 
## Residual        0.507        18          0.028                    
## Total           0.592        19                                   
## ------------------------------------------------------------------
## 
##                                  Parameter Estimates                                  
## -------------------------------------------------------------------------------------
##       model     Beta    Std. Error    Std. Beta      t       Sig      lower    upper 
## -------------------------------------------------------------------------------------
## (Intercept)    2.349         0.743                 3.162    0.005     0.788    3.910 
##         MBC    0.264         0.152        0.378    1.735    0.100    -0.056    0.584 
## -------------------------------------------------------------------------------------
result.stepre.ars.RF<-data.frame(test$predictors,test$adjr)
colnames(result.stepre.ars.RF)<-c("prd","adjr")
#View(result.stepre.ars.RF)
result.stepre.ars.RF$prd2<-c("MBC")
result.stepre.ars.RF$adjp<-c("0.100")
result.stepre.ars.RF$adjsig<-c(".")
result.stepre.ars.RF$adjr2<-c(0.1432)
library(ggplot2)
p4<-ggplot()+geom_col(data=result.stepre.ars.RF,aes(x=prd2,y=adjr),fill="#fdae61",color=NA,width = 0.5)+
  labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
  theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
                   axis.title = element_text(size = 15,colour = "black"))+
  geom_text(data = result.stepre.ars.RF,aes(x=prd2,y=adjr,label=adjsig),nudge_y = 0.05,size=7)
p4

#SARSm~
data.step<- dataRF.log[,c(5:17,19)]#dataCON[,c(3:15,17)]
library(olsrr)
model<-SARSmbc~.
fit<-lm(formula=model,data = data.step)
ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method    
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. Total_C 
## 2. Total_N 
## 3. TC_TN 
## 4. WSOC 
## 5. NH4 
## 6. NO3 
## 7. availP 
## 8. availS 
## 9. pH 
## 10. EC 
## 11. MBC 
## 12. MBN 
## 13. MBC_MBN 
## 
## We are selecting variables based on p value...
## 
## 
## Forward Selection: Step 1 
## 
## - MBC 
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.748       RMSE                0.039 
## R-Squared               0.560       Coef. Var          15.306 
## Adj. R-Squared          0.536       MSE                 0.002 
## Pred R-Squared          0.463       MAE                 0.033 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression      0.035         1          0.035    22.916     1e-04 
## Residual        0.028        18          0.002                     
## Total           0.063        19                                    
## -------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     1.083         0.173                  6.257    0.000     0.719     1.446 
##         MBC    -0.170         0.035       -0.748    -4.787    0.000    -0.244    -0.095 
## ----------------------------------------------------------------------------------------
## 
## 
## 
## No more variables to be added.
## 
## Variables Entered: 
## 
## + MBC 
## 
## 
## Final Model Output 
## ------------------
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.748       RMSE                0.039 
## R-Squared               0.560       Coef. Var          15.306 
## Adj. R-Squared          0.536       MSE                 0.002 
## Pred R-Squared          0.463       MAE                 0.033 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression      0.035         1          0.035    22.916     1e-04 
## Residual        0.028        18          0.002                     
## Total           0.063        19                                    
## -------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     1.083         0.173                  6.257    0.000     0.719     1.446 
##         MBC    -0.170         0.035       -0.748    -4.787    0.000    -0.244    -0.095 
## ----------------------------------------------------------------------------------------
## 
##                            Selection Summary                             
## ------------------------------------------------------------------------
##         Variable                  Adj.                                      
## Step    Entered     R-Square    R-Square     C(p)       AIC        RMSE     
## ------------------------------------------------------------------------
##    1    MBC           0.5601      0.5356    7.0863    -69.0077    0.0391    
## ------------------------------------------------------------------------
test<-ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method    
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. Total_C 
## 2. Total_N 
## 3. TC_TN 
## 4. WSOC 
## 5. NH4 
## 6. NO3 
## 7. availP 
## 8. availS 
## 9. pH 
## 10. EC 
## 11. MBC 
## 12. MBN 
## 13. MBC_MBN 
## 
## We are selecting variables based on p value...
## 
## 
## Forward Selection: Step 1 
## 
## - MBC 
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.748       RMSE                0.039 
## R-Squared               0.560       Coef. Var          15.306 
## Adj. R-Squared          0.536       MSE                 0.002 
## Pred R-Squared          0.463       MAE                 0.033 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression      0.035         1          0.035    22.916     1e-04 
## Residual        0.028        18          0.002                     
## Total           0.063        19                                    
## -------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     1.083         0.173                  6.257    0.000     0.719     1.446 
##         MBC    -0.170         0.035       -0.748    -4.787    0.000    -0.244    -0.095 
## ----------------------------------------------------------------------------------------
## 
## 
## 
## No more variables to be added.
## 
## Variables Entered: 
## 
## + MBC 
## 
## 
## Final Model Output 
## ------------------
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.748       RMSE                0.039 
## R-Squared               0.560       Coef. Var          15.306 
## Adj. R-Squared          0.536       MSE                 0.002 
## Pred R-Squared          0.463       MAE                 0.033 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression      0.035         1          0.035    22.916     1e-04 
## Residual        0.028        18          0.002                     
## Total           0.063        19                                    
## -------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     1.083         0.173                  6.257    0.000     0.719     1.446 
##         MBC    -0.170         0.035       -0.748    -4.787    0.000    -0.244    -0.095 
## ----------------------------------------------------------------------------------------
result.stepre.sars.RF<-data.frame(test$predictors,test$adjr)
colnames(result.stepre.sars.RF)<-c("prd","adjr")
#View(result.stepre.sars.RF)
result.stepre.sars.RF$prd2<-c("MBC")
result.stepre.sars.RF$adjp<-c("0.0001")
result.stepre.sars.RF$adjsig<-c("***")
result.stepre.sars.RF$adjr2<-c(0.5356)
library(ggplot2)
p5<-ggplot()+geom_col(data=result.stepre.sars.RF,aes(x=prd2,y=adjr),fill="#fdae61",color=NA,width = 0.5)+
  labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
  theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
                   axis.title = element_text(size = 15,colour = "black"))+
  geom_text(data = result.stepre.sars.RF,aes(x=prd2,y=adjr,label=adjsig),nudge_y = 0.05,size=7)
p5

library(ggpubr)
ggarrange(p1,p2,ncol = 1,nrow = 2,widths=c(1,1),heights=c(1,1),labels = c("(b)","(c)"))

#ggsave("plot2/Fig.2_stepwiseControls.pdf",plot = last_plot(),units = "in",height = 12,width = 6,dpi = 600)

result.stepre.ars.RF$enzyme<-"ARS-RF"
result.stepre.sars.RF$enzyme<-"SARS-RF"
result.stepre.RF<-rbind(result.stepre.ars.RF,result.stepre.sars.RF)
result.stepre.RF$forest<-"RF"
#View(result.stepre.RF)

library(ggplot2)
p6<-ggplot()+geom_col(data=result.stepre.RF,aes(x=prd2,y=adjr),fill="#fdae61",color=NA,width = 0.5)+
  labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
  theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
                   axis.title = element_text(size = 15,colour = "black"),
                   #axis.text.x = element_text(angle = 15,hjust = 1),
                   strip.text.x = element_text(size = 15,colour = "black"))+
  geom_text(data = result.stepre.RF,aes(x=prd2,y=adjr,label=adjsig),nudge_y = 0.05,size=7)+
  facet_grid(~enzyme,scales = "free")
p6

####Data DF####
dataDF.log<-subset(mydata.log,Forest=="DF")

####stepwise regression analysis####
#ARS~
data.step<- dataDF.log[,5:18]
library(olsrr)
model<-ARS~.
fit<-lm(formula=model,data = data.step)
ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method    
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. Total_C 
## 2. Total_N 
## 3. TC_TN 
## 4. WSOC 
## 5. NH4 
## 6. NO3 
## 7. availP 
## 8. availS 
## 9. pH 
## 10. EC 
## 11. MBC 
## 12. MBN 
## 13. MBC_MBN 
## 
## We are selecting variables based on p value...
## 
## 
## Forward Selection: Step 1 
## 
## - availP 
## 
##                         Model Summary                         
## -------------------------------------------------------------
## R                       0.811       RMSE               0.230 
## R-Squared               0.658       Coef. Var          6.121 
## Adj. R-Squared          0.639       MSE                0.053 
## Pred R-Squared          0.590       MAE                0.168 
## -------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression      1.837         1          1.837    34.581    0.0000 
## Residual        0.956        18          0.053                     
## Total           2.793        19                                    
## -------------------------------------------------------------------
## 
##                                  Parameter Estimates                                  
## -------------------------------------------------------------------------------------
##       model     Beta    Std. Error    Std. Beta      t        Sig     lower    upper 
## -------------------------------------------------------------------------------------
## (Intercept)    3.385         0.083                 40.902    0.000    3.211    3.559 
##      availP    0.123         0.021        0.811     5.881    0.000    0.079    0.166 
## -------------------------------------------------------------------------------------
## 
## 
## 
## Forward Selection: Step 2 
## 
## - NH4 
## 
##                         Model Summary                         
## -------------------------------------------------------------
## R                       0.850       RMSE               0.214 
## R-Squared               0.722       Coef. Var          5.676 
## Adj. R-Squared          0.689       MSE                0.046 
## Pred R-Squared          0.640       MAE                0.162 
## -------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression      2.017         2          1.008    22.079    0.0000 
## Residual        0.776        17          0.046                     
## Total           2.793        19                                    
## -------------------------------------------------------------------
## 
##                                  Parameter Estimates                                  
## -------------------------------------------------------------------------------------
##       model     Beta    Std. Error    Std. Beta      t       Sig      lower    upper 
## -------------------------------------------------------------------------------------
## (Intercept)    2.773         0.318                 8.721    0.000     2.102    3.443 
##      availP    0.121         0.019        0.802    6.268    0.000     0.080    0.162 
##         NH4    0.726         0.366        0.254    1.984    0.064    -0.046    1.498 
## -------------------------------------------------------------------------------------
## 
## 
## 
## No more variables to be added.
## 
## Variables Entered: 
## 
## + availP 
## + NH4 
## 
## 
## Final Model Output 
## ------------------
## 
##                         Model Summary                         
## -------------------------------------------------------------
## R                       0.850       RMSE               0.214 
## R-Squared               0.722       Coef. Var          5.676 
## Adj. R-Squared          0.689       MSE                0.046 
## Pred R-Squared          0.640       MAE                0.162 
## -------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression      2.017         2          1.008    22.079    0.0000 
## Residual        0.776        17          0.046                     
## Total           2.793        19                                    
## -------------------------------------------------------------------
## 
##                                  Parameter Estimates                                  
## -------------------------------------------------------------------------------------
##       model     Beta    Std. Error    Std. Beta      t       Sig      lower    upper 
## -------------------------------------------------------------------------------------
## (Intercept)    2.773         0.318                 8.721    0.000     2.102    3.443 
##      availP    0.121         0.019        0.802    6.268    0.000     0.080    0.162 
##         NH4    0.726         0.366        0.254    1.984    0.064    -0.046    1.498 
## -------------------------------------------------------------------------------------
## 
##                            Selection Summary                             
## ------------------------------------------------------------------------
##         Variable                  Adj.                                      
## Step    Entered     R-Square    R-Square     C(p)        AIC       RMSE     
## ------------------------------------------------------------------------
##    1    availP        0.6577      0.6387    -2.5829     1.9478    0.2305    
##    2    NH4           0.7220      0.6893    -3.1055    -0.2176    0.2137    
## ------------------------------------------------------------------------
test<-ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method    
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. Total_C 
## 2. Total_N 
## 3. TC_TN 
## 4. WSOC 
## 5. NH4 
## 6. NO3 
## 7. availP 
## 8. availS 
## 9. pH 
## 10. EC 
## 11. MBC 
## 12. MBN 
## 13. MBC_MBN 
## 
## We are selecting variables based on p value...
## 
## 
## Forward Selection: Step 1 
## 
## - availP 
## 
##                         Model Summary                         
## -------------------------------------------------------------
## R                       0.811       RMSE               0.230 
## R-Squared               0.658       Coef. Var          6.121 
## Adj. R-Squared          0.639       MSE                0.053 
## Pred R-Squared          0.590       MAE                0.168 
## -------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression      1.837         1          1.837    34.581    0.0000 
## Residual        0.956        18          0.053                     
## Total           2.793        19                                    
## -------------------------------------------------------------------
## 
##                                  Parameter Estimates                                  
## -------------------------------------------------------------------------------------
##       model     Beta    Std. Error    Std. Beta      t        Sig     lower    upper 
## -------------------------------------------------------------------------------------
## (Intercept)    3.385         0.083                 40.902    0.000    3.211    3.559 
##      availP    0.123         0.021        0.811     5.881    0.000    0.079    0.166 
## -------------------------------------------------------------------------------------
## 
## 
## 
## Forward Selection: Step 2 
## 
## - NH4 
## 
##                         Model Summary                         
## -------------------------------------------------------------
## R                       0.850       RMSE               0.214 
## R-Squared               0.722       Coef. Var          5.676 
## Adj. R-Squared          0.689       MSE                0.046 
## Pred R-Squared          0.640       MAE                0.162 
## -------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression      2.017         2          1.008    22.079    0.0000 
## Residual        0.776        17          0.046                     
## Total           2.793        19                                    
## -------------------------------------------------------------------
## 
##                                  Parameter Estimates                                  
## -------------------------------------------------------------------------------------
##       model     Beta    Std. Error    Std. Beta      t       Sig      lower    upper 
## -------------------------------------------------------------------------------------
## (Intercept)    2.773         0.318                 8.721    0.000     2.102    3.443 
##      availP    0.121         0.019        0.802    6.268    0.000     0.080    0.162 
##         NH4    0.726         0.366        0.254    1.984    0.064    -0.046    1.498 
## -------------------------------------------------------------------------------------
## 
## 
## 
## No more variables to be added.
## 
## Variables Entered: 
## 
## + availP 
## + NH4 
## 
## 
## Final Model Output 
## ------------------
## 
##                         Model Summary                         
## -------------------------------------------------------------
## R                       0.850       RMSE               0.214 
## R-Squared               0.722       Coef. Var          5.676 
## Adj. R-Squared          0.689       MSE                0.046 
## Pred R-Squared          0.640       MAE                0.162 
## -------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression      2.017         2          1.008    22.079    0.0000 
## Residual        0.776        17          0.046                     
## Total           2.793        19                                    
## -------------------------------------------------------------------
## 
##                                  Parameter Estimates                                  
## -------------------------------------------------------------------------------------
##       model     Beta    Std. Error    Std. Beta      t       Sig      lower    upper 
## -------------------------------------------------------------------------------------
## (Intercept)    2.773         0.318                 8.721    0.000     2.102    3.443 
##      availP    0.121         0.019        0.802    6.268    0.000     0.080    0.162 
##         NH4    0.726         0.366        0.254    1.984    0.064    -0.046    1.498 
## -------------------------------------------------------------------------------------
result.stepre.ars.DF<-data.frame(test$predictors,test$adjr)
colnames(result.stepre.ars.DF)<-c("prd","adjr")
#View(result.stepre.ars.DF)
result.stepre.ars.DF$prd2<-c("extactP","NH4+")
result.stepre.ars.DF$adjp<-c("0.0001","0.064")
result.stepre.ars.DF$adjsig<-c("***",".")
result.stepre.ars.DF$adjr2<-c(0.6387,0.0506)

library(ggplot2)
p7<-ggplot()+geom_col(data=result.stepre.ars.DF,aes(x=prd2,y=adjr),fill="#fdae61",color=NA,width = 0.5)+
  labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
  theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
                   axis.title = element_text(size = 15,colour = "black"))+
  geom_text(data = result.stepre.ars.DF,aes(x=prd2,y=adjr,label=adjsig),nudge_y = 0.05,size=7)
p7

#SARSm~
data.step<- dataDF.log[,c(5:17,19)]#dataCON[,c(3:15,17)]
library(olsrr)
model<-SARSmbc~.
fit<-lm(formula=model,data = data.step)
ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method    
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. Total_C 
## 2. Total_N 
## 3. TC_TN 
## 4. WSOC 
## 5. NH4 
## 6. NO3 
## 7. availP 
## 8. availS 
## 9. pH 
## 10. EC 
## 11. MBC 
## 12. MBN 
## 13. MBC_MBN 
## 
## We are selecting variables based on p value...
## 
## 
## Forward Selection: Step 1 
## 
## - availS 
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.594       RMSE                0.058 
## R-Squared               0.352       Coef. Var          25.110 
## Adj. R-Squared          0.316       MSE                 0.003 
## Pred R-Squared          0.210       MAE                 0.047 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.033         1          0.033    9.793    0.0058 
## Residual        0.061        18          0.003                    
## Total           0.095        19                                   
## ------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     0.664         0.139                  4.794    0.000     0.373     0.955 
##      availS    -0.131         0.042       -0.594    -3.129    0.006    -0.219    -0.043 
## ----------------------------------------------------------------------------------------
## 
## 
## 
## Forward Selection: Step 2 
## 
## - MBC 
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.694       RMSE                0.054 
## R-Squared               0.482       Coef. Var          23.103 
## Adj. R-Squared          0.421       MSE                 0.003 
## Pred R-Squared          0.271       MAE                 0.041 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.046         2          0.023    7.916    0.0037 
## Residual        0.049        17          0.003                    
## Total           0.095        19                                   
## ------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     1.403         0.380                  3.694    0.002     0.602     2.204 
##      availS    -0.153         0.040       -0.695    -3.834    0.001    -0.238    -0.069 
##         MBC    -0.129         0.063       -0.374    -2.065    0.055    -0.262     0.003 
## ----------------------------------------------------------------------------------------
## 
## 
## 
## Forward Selection: Step 3 
## 
## - availP 
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.816       RMSE                0.044 
## R-Squared               0.666       Coef. Var          19.134 
## Adj. R-Squared          0.603       MSE                 0.002 
## Pred R-Squared          0.454       MAE                 0.031 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression      0.063         3          0.021    10.621     4e-04 
## Residual        0.032        16          0.002                     
## Total           0.095        19                                    
## -------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     1.465         0.315                  4.649    0.000     0.797     2.133 
##      availS    -0.046         0.049       -0.207    -0.929    0.367    -0.150     0.059 
##         MBC    -0.223         0.061       -0.646    -3.673    0.002    -0.352    -0.094 
##      availP     0.021         0.007        0.754     2.964    0.009     0.006     0.036 
## ----------------------------------------------------------------------------------------
## 
## 
## 
## Forward Selection: Step 4 
## 
## - NH4 
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.851       RMSE                0.042 
## R-Squared               0.724       Coef. Var          17.946 
## Adj. R-Squared          0.651       MSE                 0.002 
## Pred R-Squared          0.499       MAE                 0.030 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.069         4          0.017    9.853     4e-04 
## Residual        0.026        15          0.002                    
## Total           0.095        19                                   
## ------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     1.343         0.303                  4.429    0.000     0.697     1.990 
##      availS    -0.018         0.049       -0.081    -0.369    0.717    -0.122     0.086 
##         MBC    -0.242         0.058       -0.700    -4.171    0.001    -0.365    -0.118 
##      availP     0.024         0.007        0.867     3.513    0.003     0.009     0.039 
##         NH4     0.135         0.076        0.257     1.786    0.094    -0.026     0.297 
## ----------------------------------------------------------------------------------------
## 
## 
## 
## No more variables to be added.
## 
## Variables Entered: 
## 
## + availS 
## + MBC 
## + availP 
## + NH4 
## 
## 
## Final Model Output 
## ------------------
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.851       RMSE                0.042 
## R-Squared               0.724       Coef. Var          17.946 
## Adj. R-Squared          0.651       MSE                 0.002 
## Pred R-Squared          0.499       MAE                 0.030 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.069         4          0.017    9.853     4e-04 
## Residual        0.026        15          0.002                    
## Total           0.095        19                                   
## ------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     1.343         0.303                  4.429    0.000     0.697     1.990 
##      availS    -0.018         0.049       -0.081    -0.369    0.717    -0.122     0.086 
##         MBC    -0.242         0.058       -0.700    -4.171    0.001    -0.365    -0.118 
##      availP     0.024         0.007        0.867     3.513    0.003     0.009     0.039 
##         NH4     0.135         0.076        0.257     1.786    0.094    -0.026     0.297 
## ----------------------------------------------------------------------------------------
## 
##                            Selection Summary                             
## ------------------------------------------------------------------------
##         Variable                  Adj.                                      
## Step    Entered     R-Square    R-Square     C(p)       AIC        RMSE     
## ------------------------------------------------------------------------
##    1    availS        0.3523      0.3164    8.7795    -52.9700    0.0584    
##    2    MBC           0.4822      0.4213    5.8109    -55.4457    0.0537    
##    3    availP        0.6657      0.6030    0.7898    -62.1974    0.0445    
##    4    NH4           0.7243      0.6508    0.5477    -64.0521    0.0417    
## ------------------------------------------------------------------------
test<-ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method    
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. Total_C 
## 2. Total_N 
## 3. TC_TN 
## 4. WSOC 
## 5. NH4 
## 6. NO3 
## 7. availP 
## 8. availS 
## 9. pH 
## 10. EC 
## 11. MBC 
## 12. MBN 
## 13. MBC_MBN 
## 
## We are selecting variables based on p value...
## 
## 
## Forward Selection: Step 1 
## 
## - availS 
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.594       RMSE                0.058 
## R-Squared               0.352       Coef. Var          25.110 
## Adj. R-Squared          0.316       MSE                 0.003 
## Pred R-Squared          0.210       MAE                 0.047 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.033         1          0.033    9.793    0.0058 
## Residual        0.061        18          0.003                    
## Total           0.095        19                                   
## ------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     0.664         0.139                  4.794    0.000     0.373     0.955 
##      availS    -0.131         0.042       -0.594    -3.129    0.006    -0.219    -0.043 
## ----------------------------------------------------------------------------------------
## 
## 
## 
## Forward Selection: Step 2 
## 
## - MBC 
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.694       RMSE                0.054 
## R-Squared               0.482       Coef. Var          23.103 
## Adj. R-Squared          0.421       MSE                 0.003 
## Pred R-Squared          0.271       MAE                 0.041 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.046         2          0.023    7.916    0.0037 
## Residual        0.049        17          0.003                    
## Total           0.095        19                                   
## ------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     1.403         0.380                  3.694    0.002     0.602     2.204 
##      availS    -0.153         0.040       -0.695    -3.834    0.001    -0.238    -0.069 
##         MBC    -0.129         0.063       -0.374    -2.065    0.055    -0.262     0.003 
## ----------------------------------------------------------------------------------------
## 
## 
## 
## Forward Selection: Step 3 
## 
## - availP 
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.816       RMSE                0.044 
## R-Squared               0.666       Coef. Var          19.134 
## Adj. R-Squared          0.603       MSE                 0.002 
## Pred R-Squared          0.454       MAE                 0.031 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression      0.063         3          0.021    10.621     4e-04 
## Residual        0.032        16          0.002                     
## Total           0.095        19                                    
## -------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     1.465         0.315                  4.649    0.000     0.797     2.133 
##      availS    -0.046         0.049       -0.207    -0.929    0.367    -0.150     0.059 
##         MBC    -0.223         0.061       -0.646    -3.673    0.002    -0.352    -0.094 
##      availP     0.021         0.007        0.754     2.964    0.009     0.006     0.036 
## ----------------------------------------------------------------------------------------
## 
## 
## 
## Forward Selection: Step 4 
## 
## - NH4 
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.851       RMSE                0.042 
## R-Squared               0.724       Coef. Var          17.946 
## Adj. R-Squared          0.651       MSE                 0.002 
## Pred R-Squared          0.499       MAE                 0.030 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.069         4          0.017    9.853     4e-04 
## Residual        0.026        15          0.002                    
## Total           0.095        19                                   
## ------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     1.343         0.303                  4.429    0.000     0.697     1.990 
##      availS    -0.018         0.049       -0.081    -0.369    0.717    -0.122     0.086 
##         MBC    -0.242         0.058       -0.700    -4.171    0.001    -0.365    -0.118 
##      availP     0.024         0.007        0.867     3.513    0.003     0.009     0.039 
##         NH4     0.135         0.076        0.257     1.786    0.094    -0.026     0.297 
## ----------------------------------------------------------------------------------------
## 
## 
## 
## No more variables to be added.
## 
## Variables Entered: 
## 
## + availS 
## + MBC 
## + availP 
## + NH4 
## 
## 
## Final Model Output 
## ------------------
## 
##                         Model Summary                          
## --------------------------------------------------------------
## R                       0.851       RMSE                0.042 
## R-Squared               0.724       Coef. Var          17.946 
## Adj. R-Squared          0.651       MSE                 0.002 
## Pred R-Squared          0.499       MAE                 0.030 
## --------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.069         4          0.017    9.853     4e-04 
## Residual        0.026        15          0.002                    
## Total           0.095        19                                   
## ------------------------------------------------------------------
## 
##                                   Parameter Estimates                                    
## ----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
## ----------------------------------------------------------------------------------------
## (Intercept)     1.343         0.303                  4.429    0.000     0.697     1.990 
##      availS    -0.018         0.049       -0.081    -0.369    0.717    -0.122     0.086 
##         MBC    -0.242         0.058       -0.700    -4.171    0.001    -0.365    -0.118 
##      availP     0.024         0.007        0.867     3.513    0.003     0.009     0.039 
##         NH4     0.135         0.076        0.257     1.786    0.094    -0.026     0.297 
## ----------------------------------------------------------------------------------------
result.stepre.sars.DF<-data.frame(test$predictors,test$adjr)
colnames(result.stepre.sars.DF)<-c("prd","adjr")
#View(result.stepre.sars.DF)
result.stepre.sars.DF$prd2<-c("extractS","MBC","extractP","NH4")
result.stepre.sars.DF$adjp<-c("0.717","0.001","0.003","0.094")
result.stepre.sars.DF$adjsig<-c(".","***","**",".")
result.stepre.sars.DF$adjr2<-c(0.3164,0.1049,0.1817,0.0478)
result.stepre.sars.DF<-result.stepre.sars.DF[result.stepre.sars.DF$adjp<0.05,]

library(ggplot2)
p8<-ggplot()+geom_col(data=result.stepre.sars.DF,aes(x=prd2,y=adjr),fill="#fdae61",color=NA,width = 0.5)+
  labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
  theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
                   axis.title = element_text(size = 15,colour = "black"))+
  geom_text(data = result.stepre.sars.DF,aes(x=prd2,y=adjr,label=adjsig),nudge_y = 0.05,size=7)
p8

library(ggpubr)
ggarrange(p1,p2,ncol = 1,nrow = 2,widths=c(1,1),heights=c(1,1),labels = c("(b)","(c)"))

#ggsave("plot2/Fig.2_stepwiseControls.pdf",plot = last_plot(),units = "in",height = 12,width = 6,dpi = 600)

result.stepre.ars.DF$enzyme<-"ARS-DF"
result.stepre.sars.DF$enzyme<-"SARS-DF"
result.stepre.DF<-rbind(result.stepre.ars.DF,result.stepre.sars.DF)
result.stepre.DF$forest<-"DF"
#View(result.stepre.DF)

library(ggplot2)
p9<-ggplot()+geom_col(data=result.stepre.DF,aes(x=prd2,y=adjr),fill="#fdae61",color=NA,width = 0.5)+
  labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
  theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
                   axis.title = element_text(size = 15,colour = "black"),
                   #axis.text.x = element_text(angle = 15,hjust =1),
                   strip.text.x = element_text(size = 15,colour = "black"))+
  geom_text(data = result.stepre.DF,aes(x=prd2,y=adjr,label=adjsig),nudge_y = 0.05,size=7)+
  facet_grid(~enzyme,scales = "free")
p9

library(ggpubr)
ggarrange(p3,p6,p9,ncol = 3,nrow = 1,widths=c(1,1),heights=c(1,1),labels = c("(b)","(c)","(d)"))

#ggsave("plot2/Fig.5_stepwiseControls.pdf",plot = last_plot(),units = "in",height = 4,width = 12,dpi = 600)

library(ggplot2)

p3<-ggplot(data=result.stepre.MF,aes(x=reorder(prd2, -adjr2),y=adjr2))+geom_col(fill="#fdae61",color=NA,width = 0.5)+
  labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
  theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
                   axis.title = element_text(size = 15,colour = "black"),
                   #axis.text.x = element_text(angle = 15,hjust = 1),
                   strip.text.x = element_text(size = 15,colour = "black"))+
  geom_text(data = result.stepre.MF,aes(x=prd2,y=adjr2,label=adjsig),nudge_y = 0.05,size=7)+
  geom_text(aes(label=round(adjr2,2)),vjust=2,color="black",size=3)+
  facet_grid(~enzyme,scales = "free")
p3

p6<-ggplot(data=result.stepre.RF,aes(x=reorder(prd2, -adjr2),y=adjr2))+geom_col(fill="#fdae61",color=NA,width = 0.5)+
  labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
  theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
                   axis.title = element_text(size = 15,colour = "black"),
                   #axis.text.x = element_text(angle = 15,hjust = 1),
                   strip.text.x = element_text(size = 15,colour = "black"))+
  geom_text(data = result.stepre.RF,aes(x=prd2,y=adjr2,label=adjsig),nudge_y = 0.05,size=7)+
  geom_text(aes(label=round(adjr2,2)),vjust=2,color="black",size=3)+
  facet_grid(~enzyme,scales = "free")
p6

p9<-ggplot(data=result.stepre.DF,aes(x=reorder(prd2, -adjr2),y=adjr2))+geom_col(fill="#fdae61",color=NA,width = 0.5)+
  labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
  theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
                   axis.title = element_text(size = 15,colour = "black"),
                   #axis.text.x = element_text(angle = 15,hjust =1),
                   strip.text.x = element_text(size = 15,colour = "black"))+
  geom_text(data = result.stepre.DF,aes(x=prd2,y=adjr2,label=adjsig),nudge_y = 0.05,size=7)+
  geom_text(aes(label=round(adjr2,2)),vjust=2,color="black",size=3)+
  facet_grid(~enzyme,scales = "free")
p9

library(ggpubr)
ggarrange(p3,p6,p9,ncol = 3,nrow = 1,widths=c(1,1),heights=c(1,1),labels = c("(b)","(c)","(d)"))

#ggsave("plot2/Fig.5_stepwiseForest2.pdf",plot = last_plot(),units = "in",height = 3,width = 12,dpi = 600)
#Fig. S1
#SARSmbc#
g9<-ggboxplot(mydata.log,x='Treatment',y='SARSmbc',fill = 'Treatment',facet.by = 'Forest',
              add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Specific ARS activity\nLog(nmol ug MBC-1 h-1)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g9

#SBG
g14<-ggboxplot(mydata,x='Treatment',y='SBG',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Specific BG activity\n(nmol ug MBC-1 h-1)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g14

#SXYL
g15<-ggboxplot(mydata,x='Treatment',y='SXYL',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Specific XYL activity\n(nmol ug MBC-1 h-1)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g15

#SNAG
g16<-ggboxplot(mydata,x='Treatment',y='SNAG',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Specific NAG activity\n(nmol ug MBC-1 h-1)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g16

#SLAP
g17<-ggboxplot(mydata,x='Treatment',y='SLAP',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Specific LAP activity\n(nmol ug MBC-1 h-1)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g17

#SAP
g18<-ggboxplot(mydata.log,x='Treatment',y='SACP',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Specific AP activity\nLog(nmol ug MBC-1 h-1)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g18

ggarrange(g9,g14,g15,g16,g17,g18,ncol = 2,nrow = 3,widths=c(1,1),heights=c(1,1),
          common.legend = TRUE,labels = c("(a)","(b)","(c)","(d)","(e)","(f)"))

#ggsave("plot2/Fig.3_SpecificEnzymeActivity.pdf",plot = last_plot(),units = "in",height = 15,width = 14,dpi = 600)
#Fig. S2
#SARSt
g13<-ggboxplot(mydata,x='Treatment',y='SARSt',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Specific ARS activity\n(nmol ug TC-1 h-1)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g13

#SBGt
mydata$SBGt<-mydata$BG/mydata$Total_C
g34<-ggboxplot(mydata,x='Treatment',y='SBGt',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Specific BG activity\n(nmol ug TC-1 h-1)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g34

#SXYL
mydata$SXYLt<-mydata$XYL/mydata$Total_C
g35<-ggboxplot(mydata,x='Treatment',y='SXYLt',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Specific XYL activity\n(nmol ug TC-1 h-1)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g35

#SNAG
mydata$SNAGt<-mydata$NAG/mydata$Total_C
mydata.log$SNAGt<-log(mydata$SNAGt+1)
g36<-ggboxplot(mydata.log,x='Treatment',y='SNAGt',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Specific NAG activity\nLog(nmol ug TC-1 h-1)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g36

#SLAP
mydata$SLAPt<-mydata$LAP/mydata$Total_C
mydata.log$SLAPt<-log(mydata$SLAPt+1)
g37<-ggboxplot(mydata.log,x='Treatment',y='SLAPt',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Specific LAP activity\nLog(nmol ug TC-1 h-1)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g37

#SAP
mydata$SACPt<-mydata$ACP/mydata$Total_C
mydata.log$SACPt<-log(mydata$SACPt+1)
g38<-ggboxplot(mydata.log,x='Treatment',y='SACPt',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Specific AP activity\nLog(nmol ug TC-1 h-1)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g38

ggarrange(g13,g34,g35,g36,g37,g38,ncol = 2,nrow = 3,widths=c(1,1),heights=c(1,1),
          common.legend = TRUE,labels = c("(a)","(b)","(c)","(d)","(e)","(f)"))

#ggsave("plot2/Fig.S2_SpecificEnzymeActivity2.pdf",plot = last_plot(),units = "in",height = 15,width = 14,dpi = 600)
#Fig. S3
#logS_C#
g10<-ggboxplot(mydata,x='Treatment',y='logS_C',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Log(ARS)/Log(BG+XYL)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g10

#logS_N#
g11<-ggboxplot(mydata,x='Treatment',y='logS_N',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Log(ARS)/Log(NAG+LAP)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g11

#logS_P#
g12<-ggboxplot(mydata,x='Treatment',y='logS_P',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Log(ARS)/Log(AP)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g12

#logC_N
mydata$logC_N<-log(mydata$BG+mydata$XYL)/log(mydata$NAG+mydata$LAP)
g31<-ggboxplot(mydata,x='Treatment',y='logC_N',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Log(BG+XYL)/Log(NAG+LAP)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g31

#logC_P
mydata$logC_P<-log(mydata$BG+mydata$XYL)/log(mydata$ACP)
g32<-ggboxplot(mydata,x='Treatment',y='logC_P',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Log(BG+XYL)/Log(AP)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g32

#logN_P
mydata$logN_P<-log(mydata$NAG+mydata$LAP)/log(mydata$ACP)
g33<-ggboxplot(mydata,x='Treatment',y='logN_P',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Log(NAG+LAP)/Log(AP)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g33

ggarrange(g10,g11,g12,g31,g32,g33,ncol = 2,nrow = 3,widths=c(1,1),heights=c(1,1),
          common.legend = TRUE,labels = c("(a)","(b)","(c)","(d)","(e)","(f)"))

#ggsave("plot2/Fig.4_enzymeActivity5.pdf",plot = last_plot(),units = "in",height = 15,width = 14,dpi = 600)
#Fig. S4
#TC
g19<-ggboxplot(mydata,x='Treatment',y='Total_C',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Total C (mg/g)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g19

#TN
g20<-ggboxplot(mydata,x='Treatment',y='Total_N',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Total N (mg/g)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g20

#TC_TN
g21<-ggboxplot(mydata,x='Treatment',y='TC_TN',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Total C/Total N")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g21

#WSOC
g22<-ggboxplot(mydata,x='Treatment',y='WSOC',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("WSOC (mg/kg)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g22

#NH4
g23<-ggboxplot(mydata,x='Treatment',y='NH4',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("NH4+-N (mg/kg)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g23

#NO3
g24<-ggboxplot(mydata,x='Treatment',y='NO3',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("NO3--N (mg/kg)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g24

ggarrange(g19,g20,g21,g22,g23,g24,ncol = 2,nrow = 3,widths=c(1,1),heights=c(1,1),
          common.legend = TRUE,labels = c("(a)","(b)","(c)","(d)","(e)","(f)"))

#ggsave("plot2/Fig.S3_soilproperties.pdf",plot = last_plot(),units = "in",height = 15,width = 14,dpi = 600)
#Fig. S5
#availP
g25<-ggboxplot(mydata.log,x='Treatment',y='availP',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Extractable P Log(mg/kg)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g25

#pH
g26<-ggboxplot(mydata,x='Treatment',y='pH',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("pH")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g26

#EC
g27<-ggboxplot(mydata.log,x='Treatment',y='EC',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("EC Log(us/cm)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g27

#MBC
#mydata$MBC<-mydata$MBC*1000
g28<-ggboxplot(mydata,x='Treatment',y='MBC',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("MBC (mg/g)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g28

#MBN
#mydata$MBN<-mydata$MBN/1000
g29<-ggboxplot(mydata,x='Treatment',y='MBN',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("MBN (mg/g)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g29

#MBC_MBN
g30<-ggboxplot(mydata.log,x='Treatment',y='MBC_MBN',fill = 'Treatment',facet.by = 'Forest',
               add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
  scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
  stat_compare_means(method = "t.test",comparisons = my_comparison)+
  ylab("Log(MBC/MBN)")+
  theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
        legend.text = element_text(size = size),legend.title = element_text(size = size),
        strip.text.x = element_text(size = size))
g30

ggarrange(g25,g26,g27,g28,g29,g30,ncol = 2,nrow = 3,widths=c(1,1),heights=c(1,1),
          common.legend = TRUE,labels = c("(a)","(b)","(c)","(d)","(e)","(f)"))

#ggsave("plot2/Fig.S4_soilproperties.pdf",plot = last_plot(),units = "in",height = 15,width = 14,dpi = 600)
#Fig. S6
ggplot(cor.out.all, aes(Y,X)) + 
  geom_tile(aes(fill = r), size=1)+
  geom_text(data= cor.out.all[cor.out.all$r!=0,], aes(label = r), size = 5, font="bold")+
  scale_fill_gradient2(guide = "legend", high='#2c7bb6',mid="white", low='#d7191c')+
  theme_bw()+facet_grid(~Forest)+guides(fill = guide_legend(title = 'rho'))+
  theme(axis.title= element_blank(),
        axis.text.x=element_text(colour="black", size=15, face="bold",angle = 30, hjust=1),
        axis.text.y=element_text(colour="black", size=15, face="bold"))

#ggsave("plot2/Fig.S1_spearmanCorr.pdf",plot = last_plot(),units = "in",height = 6,width = 12,dpi = 600)

R Markdown